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1.  INTRODUCTION 


Symptoms  of  decompression  sickness  (DCS)  typically  arc  assumed  to  result  from  the 
formation  of  gas  bubbles  in  blood  or  extravascular  tissue.  Information  on  extravascular 
bubbles  has  been  severely  limited  by  the  lack  of  any  noninvasive  method  of  observing 
them. 

In  a  previous  report  (NMRI  Technical  Report  91-39.  Basic  Operation  of  a  Detector 
for  Stationary  Gas  Bubbles1),  we  discussed  a  system  constructed  at  the  Jet  Propulsion 
Laboratory  (JPL)  that  uses  ultrasound  to  detect  stationary  gas  bubbles.  We  call  the 
system  the  ‘JPL  swept  frequency  bubble  detector’,  or  simply  the  ‘bubble  detector’.  The 
bubble  detector  can,  in  principle,  enable  the  quantitative  measurement  of  the  sizes  and 
numbers  of  stationary  extravascular  bubbles.  Its  theory  of  operation  depends  on  these 
facts:  1)  gas  bubbles  are  the  softest  objects  in  tissue,  and  consequently  undergo  the 
largest  deformations  when  vibrating  in  a  sound  field,  2)  although  vibrations  are  linear  at 
sufficiently  small  deformations,  they  become  increasingly  nonlinear  at  larger 
deformations2,  and  3)  nonlinear  vibrations  contain  harmonics  and  subharmonics  -  that  is, 
the  vibrations  contain  multiple  frequency  components  even  when  the  driving  signal  is 
monotonal.  In  theory  therefore,  when  tissue  is  driven  to  vibration  by  externally  applied 
sound  of  moderate  amplitude,  any  harmonics  or  subharmonics  in  the  sound 
backscattered  from  the  tissue  should  be  attributable  to  gas  bubbles  rather  than  to  stiffer 
objects.  More  specifically,  t  ie  bubble  detector  is  designed  to  exploit  the  ‘second 
harmonic’  component  of  the  energy  backscattered  by  vibrating  bubbles  (the  component 
having  the  frequency  twice  that  of  the  driving  signal). 


A  brief  description  of  the  bubble  detector  is  offered  here.  The  reader  should  look  to 
reference  1  for  a  more  comprehensive  discussion. 

The  bubble  detector  uses  a  transmitting  pressure  transducer  to  broadcast  a 
swept-frequency  sound  pressure  field  onto  a  target  and  a  receiving  pressure  transducer  to 
detect  the  sound  backscattered  from  the  target.  It  then  filters  out  all  frequency 
components  in  the  output  signal  except  for  either  the  second  harmonic  or  the 
‘fundamental’  (the  component  having  the  same  frequency  as  the  driving  signal),  as 
desired.  The  filtered  signal  is  then  processed  to  give  two  types  of  conditioned  output 
signals.  In  the  ‘range’  mode  the  system  produces  an  amplitude  spectrum  (i.e.,  a  plot  of 
amplitude  versus  frequency)  in  which  the  independent  variable  is  the  difference  in 
frequency  between  the  transmitted  and  received  signals.  This  frequency  offset  is 
proportional  to  the  time  delay  between  sound  transmission  and  reception  of  the 
backscattered  sound  from  the  target.  The  distance  from  the  transducer  he^d  to  the 
target  is  calculated  from  the  corresponding  frequency  offset. 

In  the  ‘frequency  response’  (‘FRC’)  mode,  the  system  generates  a  spectrum  in  which 
the  independent  variable  is  the  frequency  of  the  transmitted  signal  (referred  to  hence  as 
the  ‘forcing  frequency'  or  ‘driving  frequency’)  and  the  dependent  variable  is  the 
amplitude  of  the  backscattered  signal  from  any  targets  that  are  at  a  certain  preselected 
distance  from  the  transducer  head.  This  output  shows  the  ‘frequency  response’  of  the 
system  consisting  of  both  the  bubble  and  the  bubble  detector;  that  is,  it  shows  how 
strongly  the  system  responds  to  its  input  as  a  function  of  the  input  frequency.  In  theory, 
this  information  should  enable  one  to  assay  a  bubble  population:  the  location  of  a 
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resonance  peak  in  the  spectrum  depends  on  the  bubble  diameter  and  the  peak  size  is 
proportional  to  the  number  of  bubbles  at  that  diameter.  The  transducers  supplied  by 
JPL  operate  over  a  frequency  range  of  ~l-7  Mhz,  which  encompasses  the  main 
resonance  frequencies  of  bubbles  in  H20  of  roughly  0.9-6.0  microns  diameter,  so  we  can 
reasonably  expect  to  be  able  to  identify  bubbles  in  this  size  range.’ 

The  first  technical  report  in  this  series  dealt  with  the  basic  operation  of  the  bubble 
detector.  As  it  was  being  written,  the  following  goals  had  been  achieved: 

1)  A  computer  routine  that  uses  an  approximate,  analytic  solution  to  the  differential 
equations  describing  a  vibrating  bubble  to  compute  the  amplitude  of  the  second 
harmonic  component  of  the  sound  radiated  from  a  bubble  in  a  viscous  liquid  as  a 
function  of  the  frequency  of  incident  sound. 

2)  A  protocol  for  preparing  calibration  standards  consisting  of  known  bubble 
populations  by  trapping  gas  bubbles  in  transparent  hydrogels  and  determining  the 
sizes  and  numbers  of  bubbles  in  the  gels  microscopically. 

3)  Indications  that,  through  conducting  preliminary,  semi-quantitative  experiments 
with  bubbles  in  hydrogels,  the  signal/noise  ratio  of  the  system  as  delivered  was 
unacceptably  low. 

We  report  here  on  our  continuation  of  bubble  detection  studies.  The  work  generally  has 
been  directed  toward  redesigning  the  original  system  to  remove  shortcomings  that 


In  NMRI  Technical  Report  91-39,  it  was  stated  that  the  transducers  operated  over  the  frequency  range 
0.2-5  MHz,  encompassing  the  main  resonance  frequencies  of  bubbles  of  -1-20 pm  diameter.  This  statement 
was  erroneous.  By  direct  measurement,  the  frequency  resjtonsc  of  the  b  -bide  detector  system  as  a  whole  is 
reasonably  strong  over  the  range  of  - 1-7  MHz;  the  manufacturer  of  the  transducers  states  that  the 
transducers  were  designed  for  operation  over  the  2-10  MHz  range. 
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severely  limited  its  utility,  and  determining  whether  the  output  signal-to-noise  (S/N)  ratio 
is  high  enough  to  enable  quantitative  analysis  of  the  data  for  bubbles  up  to  6  /mi 
diameter.  We  eventually  had  to  evaluate  the  S/N  ratio  using  theory  rather  than 
empiricism  because  we  have  found  bubbles  this  small  to  be  too  short-lived  in  vitro  to 
permit  their  examination  using  both  the  bubble  detector  and  a  suitable  microscope  (the 
latter  being  used  to  establish  bubble  size).  In  fact,  the  theory  discussed  in  Section  V  on 
the  stability  of  bubble  populations  in  a  closed  system  suggests  that  it  is  impractical  to 
extend  the  lifetimes  of  such  small  bubbles  to  sufficient  length  for  our  purposes  in  a 
system  of  constant  surface  tension. 

The  following  new  goals  have  been  achieved  and  will  be  discussed  in  this  report: 

1)  Replacing  the  analogue  AM  tuner,  which  was  provided  with  the  original  system 
for  demodulating  AM  signals,  with  a  commercially  available  signal  analyzer 
capable  of  digital  demodulation.  This  improves  the  system’s  S/N  ratio 
somewhat.  More  importantly,  the  signal  analyzer  has  been  interfaced  with  a  PC 
upon  which  digitized  data  can  be  stored,  which  makes  the  numerical  data 
available  to  us  for  the  first  time. 

2)  Exploring  strategies  for  preparing  suitable  calibration  standards  consisting  of 
quasi-stable  bubbles  trapped  in  polymer  gels. 

3)  Elucidating  whether  small  spherical  bubbles  can  be  stabilized  (for  the  purpose 
of  producing  calibrants)  without  invoking  complicated  mechanisms  involving 
variable  surface  tension  or  variable  permeability  at  the  bubble  surface. 
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Developing  and  exploring  a  mathematical  simulation  of  the  growth  and  shrinkage 
of  bubbles  in  a  closed  system  with  finite  surface  tension. 

4)  Developing  an  approximate,  analytic  solution  to  the  differential  equations 
describing  a  vibrating  bubble  in  an  elastic  solid  and  coding  it  as  computer 
routines.  This  approach  uses  linearized  elastic  theory  (valid  for  small 
deformations)  to  model  the  elastic  stress  associated  with  oscillations  and 
nonlinear  elastic  theory  to  describe  the  potentially  large  elastic  stress  present  at 
the  bubble  wall  at  rest. 

5)  Writing  into  computer  code  a  more  accurate,  numerical  solution  to  the 
differential  equations  describing  a  vibrating  bubble  in  a  viscous  liquid  and  in  an 
elastic  solid. 

6)  Measuring  system  noise  has  been  measured. 

7)  Estimating  the  S/N  ratio  cf  the  system  for  single  bubbles  in  the  size  range  for 
which  the  system  was  designed  by  using  the  system  noise  measurements  and  the 
predictions  of  the  mathematical  models. 

The  mathematical  simulations  mentioned  in  items  4  and  5  will  allow  computation  of 
the  amplitude  of  the  2nd  harmonic  component  of  the  sound  radiated  from  a  bubble  as  a 
function  of  the  frequency  of  incident  sound. 

The  analytic  expressions  noted  in  item  4  can  be  evaluated  quickly  enough  on  even  a 
modest  computer  to  make  them  useful  for  deducing  the  sizes  and  numbers  of  bubbles  in 
a  sample  containing  an  unknown  bubble  population  using  an  iterative  curve¬ 
fitting/parameter  estimation  approach,  as  follows: 


5 


We  will  fit  a  curve  to  the  voltage  amplitude  versus  driving  frequency  data  provided 
by  the  bubble  detector  in  its  ‘frequency  response'  mode,  with  the  values  of  bubble 
diameter  and  bubble  number  taken  as  adjustable  parameters  to  be  optimized  in  the 
curve-fitting  routine.  The  data-fitting  will  consist  of  minimizing  the  sum  of  squares 
of  error  (SSE)  using  a  Marquardt  least-squares  algorithm.3  The  SSE  is  a  measure 
of  the  goodness-of-fit  of  the  model  to  the  data.  The  independent  variable  is  the 
driving  frequency,  the  dependent  variable  is  the  voltage  amplitude,  and  for  an 
unknown  target  sample  the  parameters  to  be  estimated  are  bubble  radius  and 
bubble  number.  The  Marquardt  algorithm  is  an  iterative  routine  that 
systematically  adjusts  the  values  of  each  unknown  parameter  ‘B’  based  on  the  value 
of  <?(SSE)/dB  in  the  current  iteration. 

The  numerical  solutions  mentioned  in  item  5  are  much  more  CPU-intensive  than  the 
analytic  expressions,  but  presumably  more  accurate,  and  will  be  used  to  evaluate  the 
accuracy  of  the  analytic  solutions.  There  is  a  lack  of  experimental  verification  of  any 
theory’s  predictions  regarding  the  2nd  harmonic  component  of  bubble  oscillations. 

The  simulation  of  bubble  dynamics  (listed  as  item  3  above)  indicates  that  quasi-stable 
bubbles  of  diameters  <  100  pm  cannot  be  obtained  in  a  system  of  constant  surface 
tension.  Our  own  efforts  to  obtain  small,  quasi-stable  bubbles  (mentioned  above  in  item 
2)  produced  no  results  to  the  contrary:  we  were  unsuccessful  in  attempts  to  prepare 
suitable  calibration  standards  containing  bubbles  in  the  size  range  of  interest. 

The  mathematical  descriptions  of  bubble  behavior  from  items  4  and  5  enable  us  to 
estimate  the  output  voltage  amplitude  that  should  be  observed  for  a  6  |im  bubble  driven 
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at  its  main  resonance  frequency.  (This  is  the  strongest  output  signal  that  would  be 
observed  ^r  any  bubble  of  ^  6  um  diameter  at  any  driving  frequency.)  By  comparison 
with  the  measured  noise  floor,  we  have  estimated  the  S/N  ratio.  Ibe  result  of  this 
estimate  casts  some  doubt  on  the  practicality  of  quantitatively  studying  individual  bubbles 
in  the  0.9-6  /um  range.  The  proper  course  of  action  may  be  to  study  larger  bubbles,  thus 
increasing  the  S/N  ratio  without  improving  any  electronic  components.  This  necessitates 
the  purchase  of  pressure  transducers  designed  to  operate  over  a  lower  frequency  range 
than  the  pair  originally  supplied.  The  size  range  of  stationary  bubbles  that  are  most 
relevant  to  DCS  is  unknown,  so  there  is  no  reason  to  insist  upon  studying  bubbles  of  any 
particular  size  range  so  long  as  the  size  is  physiologically  plausible. 


II.  ADDING  DIGITAL  AM  SIGNAL  DEMODULATION  TO  THE  SYSTEM 


A.  Impetus  for  adding  digital  AM  demodulation 

As  noted  in  the  Introduction,  the  bubble  detector  has  two  operating  modes:  the 
‘range’  mode,  in  which  the  signal  is  Fourier-transformed  to  give  information  on  bubble 
location,  and  the  ‘frequency  response’,  or  FRC  mode,  in  which  the  signal  is  demodulated 
to  give  information  about  the  sizes  and  numbers  of  bubbles. 

For  the  frequency  response  mode,  the  bubble  detector  (as  delivered  by  JPL)  included 
a  section  for  analogue  demodulation  of  amplitude  modulated  (AM)  signals;  i.e.  an  AM 
tuner.  The  tuner  consisted  of  a  local  oscillator,  a  mixer,  and  a  low-pass  filter.  It  was 
noted  in  reference  1  that  the  tuner  section  was  a  significant  noise  source.  It  was  noted 
also  that  the  bubble  detector’s  output  signal  was  analog  and  that  no  means  had  been 
provided  for  digitizing  or  recording  the  output,  preventing  meaningful  data  analysis.  In 
addition,  the  output  was  not  an  amplitude  spectrum,  as  expected  from  a  system  intended 
to  function  as  a  spectrum  analyzer,  but  rather  was  the  sum  of  waveforms  of  various 
freque1  ides,  each  waveform  associated  with  one  of  the  acoustic  targets  in  front  of  the 
pressure  transducers.  This  sort  of  output  does  not  lend  itself  readily  to  objective 
analysis. 

In  order  for  the  bubble  detector  to  have  potential  as  a  research  tool,  the  problems 
noted  above  must  be  overcome.  Replacing  the  AM  tuner  with  a  digital  signal  analyzer 
capable  of  AM  signal  demodulation  can  solve  all  three  of  them.  Digital  signal 
processing  does  not  contribute  appreciably  to  noise,  commercially  available  digital 
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processors  normally  can  be  interfaced  with  computers  on  which  data  can  be  stored  and 
processed,  and  we  will  see  that  the  output  of  an  AM  demodulator  for  this  system  is  an 
amplitude  spectrum  whose  relationship  with  the  physical  system  under  study  is  easy  to 
conceptualize. 

B.  Hardware 

Demodulating  an  AM  signal  is  equivalent  to  recovering  the  amplitude  versus  time 
information;  for  this  system  the  input  frequency  is  a  linear  function  of  time,  and 
therefore  the  demodulated  signal  is  equivalent  to  a  sj.etrum  of  amplitude  versus  forcing 
frequency.  A  Hewlett-Packard  3561 A  signal  analyzer  was  chosen  for  the  task.  This 
machine  cannot  actually  demodulate  signals,  but  it  can  separate  a  signal  into  its  real  and 
imaginary  parts.  These  are  captured  and  sent  to  a  FC,  on  which  the  amplitude  is 
computed  at  each  time  point  simply  by  taking  the  square  root  of  the  sum  of  the  squares 
of  the  real  and  imaginary  parts.  The  HP  3561A  simultaneously  carries  out  a  fast  Fourier 
transformation  (FFT)  on  the  signal,  thus  performing  the  range  mode  task  at  the  same 
time  as  it  handles  the  FRC  mode  function. 

C.  jflfiwaig 

The  Appendix  contains  the  HP  Basic  program  used  to  remotely  control  the  HP  3561A 
and  the  HP  3325A  signal  generator  from  an  IBM-compatible  PC  in  which  a  Hewlett- 
Packard  Basic  Language  Processor  card  has  been  installed.  The  program  runs 
measurement  ‘cycles’,  pausing  between  cycles  for  a  user-selectable  length  of  time. 

During  each  cycle  the  program  1)  prompts  the  bubble  detector  to  make  a  series  of 
measurements,  2)  extracts  the  raw  output  from  these  measurements  from  the  signal 
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analyzer,  3)  converts  them  to  amplitude-versus-forcing  frequency  spectra  for  the  FRC 
output,  4)  averages  the  spectra  from  multiple  measurements  (i.e.,  uses  ‘signal  averaging’) 
to  minimize  white  noise,  5)  subtracts  from  the  spectrum  a  signal-averaged  "background" 
measurement  taken  at  the  beginning  of  the  experiment  when  bubbles  are  absent,  6) 
stores  the  result  to  hard  disk,  and  7)  plots  the  result  on  a  Hewlett-Packard  printer.  The 
user  chooses  an  appropriate  target  for  the  background  measurement.  For  example,  if 
bubbles  are  to  be  formed  in  an  acoustical  target  during  the  experiment,  one  chooses  as 
background  the  target  itself  before  the  introduction  of  bubbles. 

Simultaneously,  the  3561 A  also  performs  an  FFT  on  the  data  for  the  range  mode 
output.  The  result  is  transferred  to  the  PC  along  with  the  FRC  spectrum  and  is  then 
printed.  A  few  examples  of  the  final  outputs  from  both  the  FRC  and  range  modes  are 
shown  in  the  following  section. 

The  program  allows  the  user  to  select  values  of  the  following  operating  parameters: 

1)  frequency  range  of  the  frequency  sweep  (i.e.,  the  start  and  stop  frequencies) 

2)  sweep  time 

3)  voltage  amplitude  generated  by  the  3325A  signal  generator 

4)  analyzer  frequency  ‘span’  (explained  below) 

5)  ‘start  frequency’  for  the  analyzer  (explained  below) 

6)  number  of  sweeps  to  be  averaged  together  during  the  cycle  in  which  the 
background  noise  is  measured 

7)  number  of  sweeps  to  be  averaged  together  during  each  cycle  in  which  a  signal  is 
measured 
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8)  how  long  to  pause  between  cycles 
D.  On  the  selection  of  operating  parameters 

We  will  now  consider  the  analyzer’s  start  frequency,  its  frequency  span,  and  how  to 
chose  each  correctly.  The  HP  3561 A  analyzer  captures  all  signals  having  frequencies 
between  the  start  frequency  and  the  sum  of  [start  frequency  +  frequency  span].  Higher 
and  lower  frequencies  are  filtered  out.  Therefore,  the  span  is  the  bandwidth  of  the 
captured  signals.  The  user  wants  the  signal  analyzer  to  capture  the  output  signal  only 
while  the  signal  generator  is  performing  a  frequency  sweep.  The  synchrony  between 
these  two  instruments  is  realized  by  having  the  analyzer  triggered  by  the  generator  at  the 
start  of  the  sweep.  The  analyzer  also  should  stop  capturing  data  at  the  same  time  as  the 
sweep  is  completed.  In  other  words,  the  ‘time  record’  captured  by  the  analyzer  should 
be  the  same  duration  as  the  sweep  time.  The  correct  choice  of  span  ensures  this,  as  we 
discuss  now. 

The  HP  3561 A  analyzer  always  collects  210=1024  digital  samples  per  time  record. 

The  sample  rate  is  thereby  determined: 

sample  rate  =  (1024  samples)/(length  of  time  record)  [1] 

The  sample  rate  is  always  2.56  •  span.  Because  the  time  record  length  should  equal  the 
sweep  time,  we  see  by  inspection  that 

span  =  400/(sweep  time)  [2] 
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where  the  span  Is  in  Hz  and  the  sweep  time  is  in  seconds.  The  program  automatically 
selects  the  correct  span  for  a  given  sweep  time,  and  vice-versa.  One  must  manually 
select  the  start  frequency  such  that  signals  of  interest  have  frequencies  greater  than  the 
start  frequency  and  less  than  [start  frequency  +  span]. 

As  an  illustration,  the  frequency  of  the  signal  received  from  an  acoustical  target  is 
given  by 

2  (bandwidth  of  sweep,  Hz)  •  (distance,  cm) 

^signal  =  ~  PI 

(sweep  time,  sec)  •  (velocity  of  sound,  ~  1.5  •  l(r  cm/sec  in  water) 

where  the  ‘distance’  in  the  numerator  is  measured  between  the  transducer  head  and  the 
target1.  For  example,  a  0.2-second  sweep  from  1  MHz  to  7  MHz  yields  a  signal  of  1.3 
kHz  for  3,  bubble  located  3.25  cm  from  the  transducers  (the  transducer  head  is  shaped  to 
optimize  performance  for  a  distance  of  roughly  3.25  cm).  If  we  record  a  time  record  for 
a  0.2-second  sweep,  the  frequency  span  will  be  2  kHz  (as  shown  in  Equation  [2]).  We 
might  set  the  start  frequency  at  1  kHz  to  capture  all  signals  between  0.8  and  2.8  kHz, 
thus  ensuring  that  the  target  associated  with  the  1.3  kHz  signal  is  detected. 

In  the  FRC  mode  any  signals  within  the  frequency  span  are  lumped  together  in  the 
conversion  to  amplitude,  so  that  there  is  no  reso'ution  of  frequencies  for  any  signals  that 
fall  within  the  frequency  span.  Therefore,  for  the  example  given  there  is  no  frequency 
resolution  in  the  FRC  mode  for  signals  ranging  from  0.8  to  2.8  kHz,.  The  most  precise 
statement  we  could  make  about  a  bubble’s  location,  based  only  on  the  information  from 
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the  FRC  mode  output,  would  be  that  its  distance  from  the  transducer  head  is  between 
2  cm  and  7  cm. 

Although  the  range  mode  does  give  us  the  information  necessary  for  deducing  the 
distances  of  bubbles  from  the  transducers,  the  presence  of  more  than  one  bubble  will 
lead  to  ambiguity  about  which  peak  in  the  range  mode  spectrum  corresponds  to  which 
peak(s)  in  the  FRC  mode.  In  other  words,  associating  bubble  sizes  with  bubble  locations 
may  be  a  hit-or-miss  procedure. 

One  further  point:  whereas  there  is  a  large  DC  noise  component  and  also  some  noise 
at  60  Hz,  the  start  frequency  should  be  set  no  lower  than  a  few  hundred  Hz  so  that  these 
spurious  components  are  digitally  filtered  out  by  the  3561 A  before  further  signal 
processing. 


III.  THE  OUTPUT  SIGNAL 


In  this  section  examples  of  the  system’s  output  signal  will  be  shown  and  discussed  for 
various  acoustical  targets.  In  all  cases  the  amplitude  of  the  2nd  harmonic  component  of 
the  backscattered  signal  is  plotted. 

All  measurements  in  this  section  were  made  under  the  following  conditions: 
sweep  time  =  0.2  sec 

span  =  2  kHz  (therefore,  time  record  length  =  0.2  sec) 
sweep  start  fr  equency  =  2  MHz 
sweep  stop  frequency  =  10  MHz 
windowing  =  flattop 

signal  averaging:  15  time  records  averaged 
target  distance  -  3.25  cm 

The  frequency  range  is  that  over  w'hich  the  transducer  pair  was  designed  to  operate, 
not  necessarily  the  optimal  range.  The  measurements  were  made  using  signal  averaging 
over  15  frequency  sweeps.  It  was  observed  that  the  averaged  signal  changes  negligibly 
after  -5  averages,  and  therefore  after  15  averages  the  uncorrelated  noise  has  been 
minimized,  so  that  15  averages  yields  the  same  result  as  would  averaging  indefinitely. 

Figures  la  and  lb  show  measurements  made  when  the  transmitted  signal  is  reflected 
from  a  slab  of  Wall-Gone  sound  absorber.  The  upper  plots  are  the  FRC  mode  outputs 
and  the  lower  ones  are  the  range  mode  outputs  (explained  in  the  Introduction).  The 
FRC  portion  of  Figure  la  is  called  a  ‘background’  measurement  and  the  FRC  portion  of 
Figure  lb  is  a  measurement  of  [signal  +  background  -  background].  There  is  no 
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Figure  1a:  Output  signal;  Acoustical  target  is 
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background  subtraction  being  performed  on  the  range  mode  output.  Because  Figure  lb 
was  obtained  simply  by  repeating  the  background  measurement,  its  FRC  output  should 
show  nearly  zero  voltage,  as  it  does.  The  time  at  which  the  signal  was  measured  is 
recorded  in  the  header  of  each  plot. 

In  Figure  2,  we  see  similar  measurements  made  with  a  3-mm-thick  slab  of 
polyacrylamide  (PAAm)  gel  sitting  on  top  of  the  Wall-Gone.  There  are  no  bubbles  in 
the  gel. 

Under  each  FRC  plot  we  record  the  "positive  area"  between  the  x-axis  and  the 
positive  portion  of  the  curve  and  the  "negative  area"  between  the  x-axis  and  the  negative 
portion  of  the  curve,  computed  using  trapezoidal  integration.  They  are  useful  indices 
that  quantify  the  results  in  an  easily  understood  way. 

Figures  lb  and  2b  demonstrate  that  subtracting  an  appropriate  background  can  reduce 
the  noise  by  at  least  an  order  of  magnitude.  It  is  seen  that  the  residual  after  background 
subtraction  is  much  greater  with  the  PAAm  target  (Figure  2)  than  with  the  Wall-Gone. 
This  observation  has  been  reproduced  in  several  additional  measurements  with  both 
PAAm  and  Wall-Gone  (not  shown).  Obviously,  the  backscattered  signal  changes 
appreciably  from  one  measurement  cycle  to  the  next  when  replicate  measurements  are 
made  with  a  PAAm  gel.  The  reason  for  this  change  is  not  yet  clear. 

The  lumpiness  of  the  frequency  responses  seen  in  the  background  measurements 
in  Figures  1  and  2  results  from  beating  between  two  or  more  signals  of  similar 
frequencies.  To  understand  this,  consider  that  in  the  experiment  of  Figure  2  the  bubble 
detector  is  receiving  signals  reflected  from  both  the  front  and  back  surfaces  of  the 
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Figure  2a:  Output  signal;  Acoustical  target  is 
a  3  mm-thick  PAAm  gel  mounted  on 
Wall-Gone;  Background  measurement. 
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Figure  2b:  Output  signal;  Acoustical  target  is 
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PAAm  gel  simultaneously.  Because  the  signal  from  the  back  surface  has  travelled  2  •  3 
mm  =  6  mm  farther  (or  ~0.4  /isec  longer),  its  offset  frequency  should  be  -160  Hz 
higher.  The  range  mode  output  does  in  fact  show  two  distinct  frequency  components 
separated  by  160  Hz.  This  is  taken  to  mean  that  the  signal  analyzer  is  receiving  two 
signals  that  differ  in  frequency  by  160  Hz.  When  they  arc  summed  together  for  the  FRC 
display,  the  beat  phenomenon  appears  as  a  160-Hz  oscillation  in  the  amplitude.  Since 
the  time  of  one  measurement  was  0.2  second,  about  32  of  these  beats  should  appear  in 
the  FRC  plot,  as  is  the  case. 
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IV.  ATTEMPTS  TO  PREPARE  SUITABLE  CALIBRATION  STANDARDS 


For  calibrating  the  bubble  detector,  it  is  necessary  to  have  standards  that  are  subject 
to  assay  by  some  reliable  independent  method.  We  have  prepared  bubble  samples  of  the 
sort  that  we  anticipate  using  as  calibrants  and  used  them  in  some  preliminary 
measurements  intended  to  estimate  the  bubble  detector’s  S/N  ratio.  We  plan  to  use 
calibrants  consisting  of  bubbles  trapped  in  transparent  hydrogels  and  to  assay  them  using 
differential  interference  microscopy.  Our  approach  is  described  below. 

A.  Methods 

To  produce  the  gels,  aqueous  solutions  were  prepared  in  the  following  compositions: 

solution  A 

0.20  g/cm3  acrylamide  monomer 

0,01  g/cm3  N,N '  -methylene-bis-acrylamide  crosslinking  monomer 

1.0  n  1/cm3  N,N,N  ',N  '-tetramethylethylenediamine  accelerator 

solution  B 

variable  0.2-6  mg/cm3  ammonium  persulfate  initiator 

Prior  to  mixing  the  solutions,  each  was  sparged  with  nitrogen  to  remove  dissolved 
oxygen,  which  inhibits  many  vinyl  polymerizations4.  Each  solution  was  passed  through  a 
filter  having  a  0.22-jum  nominal  pore  size  to  remove  particles  that  otherwise  would 
severely  confuse  the  microscopic  examination.  Equal  amounts  of  solutions  ‘A’  and  ‘B’ 
then  were  mixed  and  the  resulting  solution  was  cast  between  either  a  pair  of  cleaned 
glass  plates  or  a  glass  microscope  slide  and  a  glass  cover  slip.  In  experiments  during  the 
period  covered  by  this  report,  the  plates  (or  the  microscope  slide  and  cover  slip)  were 
separated  by  a  150  yim-thick  plastic  spacer  designed  for  use  in  casting  electrophoresis 
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gels.  The  monomer  solution  ‘sandwich’  was  placed  under  vacuum  for  a  few  minutes 
(variable  length  of  time).  The  polymerization  was  then  allowed  to  proceed  at  ambient 
pressure  at  least  until  gelation  occurred.  Bubbles  that  formed  by  cavitation  at  the 
reduced  pressure  were  trapped  in  the  finished  gel. 

In  some  experiments  the  monomer  solution  was  chilled  to  0  °C  and  maintained  at  that 
temperature  before  and  during  sparging  with  N2,  to  increase  the  concentration  of 
dissolved  N2.  Thus,  the  overpressure  in  the  monomer  solution,  defined  as  [partial 
pressure  of  dissolved  gas]  +  [hydrostatic  pressure],  would  be  higher  under  vacuum  and 
would  persist  even  after  the  sample  had  been  returned  to  ambient  hydrostatic  pressure. 

In  other  experiments,  a  detergent  (Tween  40)  was  added  to  the  monomer  solution  to 
reduce  the  surface  tension. 

In  order  to  be  useful,  the  calibration  standards  must  be  the  same  size  when  they  are 
interrogated  ultrasonically  as  when  they  were  assayed  microscopically.  Accordingly,  the 
stability  of  these  bubbles  over  time  was  tested  by  examining  some  of  them  under  a  bright 
field  microscope  over  several  hours’  time  while  leaving  the  gel  membrane  sandwiched 
between  the  glass  plates.  Our  plan  was  to  assay  the  calibration  standards  microscopically 
and  use  the  bubble  detector  to  interrogate  them  within  24  h  of  their  production;  during 
that  time  their  dehydration  would  be  prevented  by  storage  in  an  enclosed  container 
whose  interior  was  kept  at  the  dew  point  by  the  presence  of  standing  water. 

B.  Results 

Briefly,  free  bubbles  that  initially  were  smaller  than  50  fx m  in  diameter  survived  no 
more  than  a  few  hours  before  disappearing  in  any  of  the  experiments.  Bubbles  larger 


22 


than  100 /im  survived  for  many  hours,  although  their  sizes  were  not  measured  accurately 
enough  over  time  to  allow  a  more  precise  statement  about  thtir  size  stability.  Adding 
detergent  to  the  monomer  solution  or  keeping  the  monomer  solution  chilled  during 
sparging  with  N2  dramatically  increased  the  number  of  bubbles  formed,  but  neither 
treatment  had  the  effect  of  stabilizing  the  small  bubbles,  over  several  hours.  Bubbles 
formed  spontaneously  at  ambient  pressure  in  the  chilled  monomer  solution  upon  its 
removal  from  the  ice  bath. 

The  shrinkage  and  disappearance  of  bubbles  was  easily  viewed  at  magnification  100  x . 
The  collapsing  bubbles  left  behind  spherical,  water-filled  cavities  in  the  gel  matrix  that 
were  clearly  visible  as  discontinuities.  By  staining  the  gel  with  food  coloring  we  could  be 
sure  of  which  spherical  cavities  contained  gas  and  which  contained  water  —  the  gas 
phases  transmitted  more  light  since  the  dye  was  excluded  from  them. 

In  some  experiments,  a  thin  film  of  bone  wax  was  smeared  onto  one  spot  on  a 
microscope  slide  before  the  monomer  solution  was  cast  onto  the  slide.  Bone  wax  is  very 
hydrophobic  and  was  found  to  be  an  excellent  nucleation  site,  as  expected.  Even  small 
bubbles  persisted  many  hours  on  the  bone  wax.  This  was  expected  because  the  bone  wax 
presents  an  irregular  surface  with  numerous  crevices,  and  it  is  known  that  some  bubbles 
in  crevices  can  remain  static  indefinitely5.  However,  they  are  unsuitable  for  use  as 
calibration  standards  because  no  reasonably  complete  mathematical  description  of  the 
nonlinear  dynamics  of  either  non-spherical  bubbles  or  of  bubbles  at  u  solid  boundary  is 
available. 


23 


V.  SIMULATION  OF  THE  GRO  WTH  AND  SHRINKAGE  OF  BUBBLES  IN  A 
CLOSED  SYSTEM  WITH  FINITE  SURFACE  TENSION 

A  requisite  for  calibration  standards  for  the  bubble  detector  is  that  the  bubbles  must 
remain  constant  in  size  for  enough  time  to  permit  their  examination  under  a  microscope 
and  their  ultrasonic  interrogation.  This  means  that  they  must  be  stable  over  a  period  of 
at  least  a  few  hours. 

To  give  us  an  idea  of  how  such  stable  calibrants  might  be  prepared,  a  model  was 
developed  that  predicts  bubble  diameter  changes  as  a  function  of  time  in  a  closed 
system.  Ideally,  one  would  do  this  by  solving  the  unsteady  state  continuity  equation  with 
radial  symmetry.  However,  the  continuity  equation  for  this  system  is  partial  differential 
equation  (PDE)  because  of  the  spatial,  as  well  as  temporal,  variation  in  solute 
concentration.  Solving  a  PDE  was  felt  to  be  more  time-consuming  than  this  analysis 
warranted.  Instead,  we  simplified  the  mathematics  considerably  by  assuming  that  the 
flux  of  gas  from  the  bubble  surface  is  proportional  to  the  mass  transfer  coeffic'ent 
evaluated  for  a  sphere  in  a  quiescent  liquid  at  steady-state  (see  Equations  [4]  and  [tl] 
below).  The  assumption  of  steady-state  implies  that  the  radius  of  the  sphere,  the  solute 
concentration  at  the  surface  of  the  sphere,  and  the  solute  concentration  far  from  the 
sphere,  all  are  time-invariant.  We  also  neglected  the  convective  flux  associated  with  bulk 
flow  toward  or  away  from  the  sphere  as  its  size  changes.  All  of  the  above  simplifications 
are  justified  when  the  bubble  size  is  changing  sufficient  slowly,  and  therefore  are  suitable 
for  a  study  whose  goal  was  to  identify  conditions  leading  to  time-invariant  bubble  sizes. 


24 


Thermal,  viscous,  and  compressibility  effects  were  ignored,  which  invalidates  the 
model  at  very  small  diameters  (i.e.,  below  10  nm  for  a  surface  tension  of  73  dyn/cm), 
when  implosion  is  rapid  and  internal  pressures  are  high.  This  is  inconsequential  for  our 
application:  first,  we  are  uninterested  in  knowing  precisely  when  a  fast-shrinking  bubble 
disappears;  second,  the  system  in  toto  is  unaffected  by  these  very  small  bubbles  because 
they  contain  too  little  gas  to  be  of  consequence.  The  partial  pressure  of  water  in  the  gas 
phase  is  taken  to  be  the  vapor  pressure  of  water  at  ambient  temperature.  Its 
dependency  on  either  hydrostatic  pressure  and  temperature  is  ignored,  but  again,  these 
effects  are  significant  only  for  tiny  bubbles. 

A.  Mathematical  details 

For  a  single  spherical  bubble  the  mole  flux  of  gas  away  from  its  surface  is 

J  =  k(Ciurf-C°°  )  [4] 

where  J  =  mole  flux,  mol/(cm2-sec); 

k  =  external  mass  transfer  coefficient,  cm/sec; 

C  surf  =  concentration  of  dissolved  gas  at  the  outer  bubble  surface,  mol/cm3; 

C 00  =  concentration  of  dissolved  gas  far  from  the  surface,  mol/cm3. 

We  now  make  the  following  reasonable  assumptions: 

a.  There  is  equilibrium  between  phases  across  the  phase  interface. 

b.  There  is  no  resistance  to  mass  transport  in  the  gas  phase,  so  that  the  partial 
pressure  of  gas  at  the  inner  bubble  surface  "P  turf"  simply  equals  its  bulk  value  "P  bub"  in 
the  bubble. 

c.  The  concentration  of  dissolved  gas  ^dis  is  related  to  the  partial  pressure  of  gas  with 
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which  it  would  be  in  equilibrium  P*s  according  to  a  concentration-independent  partition 
coefficient  K  (i.e.,  Cdis  =  KPdis),  which  is  a  statement  of  Henry’s  Law. 

Equation  [4]  now  becomes 

J  =  k  K  (  P  bub  -  P  00  )  [5] 

We  further  note  that  Pbub  equals  the  difference  between  the  hydrostatic  pressure  within 
the  bubble  and  the  vapor  pressure  of  H20,  and  we  account  for  the  surface  tension  effect: 

Pbub  =  9  +  2o/R  -  PH20  [6] 

where  P  =  ambient  hydrostatic  pressure,  dyn/cm2; 

rr  =  curfary*  t^ncinn  rlvn  /rm 
v  •«*’ ,  w — 

The  mass  balance  on  the  gas  inside  the  bubble  (assuming  sphericity)  gives  us 

J  =  -1/(4kR2)  dn/dt  [7] 

where  R  =  bubble  radius,  cm; 

n  =  number  of  moles  of  gas  inside  the  bubble,  mo). 

Assuming  ideal  gas,  we  have 

n  =  [p  bub(4/3)nR3]/(RGT) 


=  [  (  P  +  2o/R  -  PH2Q  )  (4/3)n  R3  ]/(RgT) 


[8] 


and 


dn/dt  =  4*/(RgT)  [  (P-P„20)R2  +  (4/3)Ro  ]  dR/dt 


[9] 


where  RG  -  gas  law  constant,  82,057  atm-cm3/mol-°K; 

T  =  absolute  temperature,  °K. 

We  can  now  combine  Equations  [5],  [6],  [7],  and  [9]  to  yield  an  ordinary  differential 
equation  for  the  rate  of  change  of  the  bubble  radius: 

-  kKRcT  (  9  -  PH20  +  2o/R  -  P°°  ) 

dR/dt  = - [10] 

[  P  *  PH2o  +  (4/3)o /R  ] 

To  compute  the  mass  transfer  coefficient  we  use  the  result  that,  for  a  sphere  in  a 
quiescent  liquid  at  steady-state  the  dimensionless  Sherwood  Number  is  2.06,  so  that 

Sh  =  2.0  =  2  R  k  /  D 

or  k  =  D/R  [11] 

where  D  is  the  diffusivity  (cm2/sec).  Substitution  into  Equation  [10]  gives 

-  DKRgT  (  P  -  Ph2o  +  2o/R  -  P  00  ) 

dR/dt  - -  [12] 

R  [  P  -  PH2o  +  (4/3)0  /R] 


Or,  when  more  than  one  bubble  is  present,  we  can  write  Equation  [12]  in  this  form  for 
bubble  i: 


27 


-  DKRgT  (  P  -  FlI20  +  2a/Rj  -  P  00  ) 

dRj/dt  = - 

Ri  [  p  -  PH2o  +  (4/3)o /RJ 


[13] 


Equation  [12]  or  [13j  is  easily  solved  analytically  for  the  special  case  of 

constant  P”,  that  is,  an  infinitely  large  system  in  which  the  bubbles  have  no  effect  on  the 

bulk  concentration  of  dissolved  gas.  For  the  finite  closed  system,  a  mass  balance  shows 

that 

dC°°/dt  =  KdP^/dt  =  -(1/VJ  E  dr^/dt  [14] 

for  i  =  1  to  N  bubbles,  where  VL  =  the  volume  of  the  condensed  phase,  in  cm3.  Then 
for  each  nubble  i,  Equation  [9]  can  be  substituted  into  Equation  [14]  to  yield 

-K  dP  ~/dt  -  4k  /(VlR0T)  E  (  [  (P-Pmt)R,2  +  4/3R,«  ]  dRj/dt  )  [15] 


This  is  integrated  to  yield 

p»  =  P0“-4.t/(3KV1.RGT)[(P-PH20)(ERj5-ER<l/)  +  2o(ERj2-ER1,/)]  [16] 

where  P0  00  =  initial  partial  pressure  of  dissolved  gas  in  the  condensed  phase,  atm. 

All  summations  are  from  i  =  1  to  N,  where  N  is  the  number  of  bubbles. 
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Therefore,  Equations  [13]  and  [16]  define  this  system;  Equation  [13]  must  be  solved 
for  each  bubble  in  turn,  simultaneously  with  Equation  [16],  We  did  this  using  a  4th- 
order  Runge-Kutta  numerical  routine. 

B.  Simulation  results 

Understanding  the  simulation  results  is  easier  if  one  keeps  in  mind  the  definition  of 
‘critical  radius’.  This  is  the  radius  of  a  bubble  that  would  be  in  both  mechanical  and 
thermodynamic  equilibrium  with  the  condensed  phase,  so  that  there  would  be  no  mass 
flux  at  its  surface  and  its  size  would  not  change  with  time.  To  derive  an  expression  for 
Rcrj„  remember  that  the  driving  force  for  mass  transport  is  (Pbub  -  P“),  or 

driving  force  =  (  9  +  2o/R  -  PH20  -  P°°  )  [17] 

A  quasi-stable  bubble  thus  can  be  defined  mathematically  as  one  for  which  the  right  side 
of  Equation  [17]  is  near  zero;  that  is,  Pbub  =  P®.  R  equals  the  critical  radius  when  the 
driving  force  is  zero,  so  we  can  solve  Equation  [12]  for  the  critical  radius: 

^crit  =  2o/(  P  +  Ph20  )  [18] 

When  R>Rcnt  the  bubble  grows  and  when  R<Rcril  it  shrinks.  In  an  unbounded  system 
where  P 00  ■  constant  the  growth  or  shrinkage  is  monotonic.  In  a  finite  system  in  which 
P®  is  a  time-dependent  variable,  a  bubble  may  alternately  grow  and  shrink  as  F® 
fluctuates  and  Rrri,  fluctuates  with  it. 
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To  further  understand  the  simulation  results,  it  should  be  noted  that  two  bubbles  of 
different  sizes  cannot  simultaneously  be  equilibrated  with  the  same  condensed  phase 
because  their  internal  pressures  differ,  and  that  a  larger  bubble  must  always  be  either 
growing  faster  or  shrinking  slower  than  a  smaller  bubble.  These  observations  lead  to  the 
conclusion  that  any  closed  system  must  ultimately  reach  one  of  two  possible  final  states: 

1)  all  of  the  bubbles  have  collapsed,  or 

2)  only  one  bubble  remains,  and  it  is  the  one  that  was  largest  initially. 

So,  we  see  that  bubbles  initially  grow  or  shrink  depending  upon  whether  their  radii  are 
greater  than  or  less  than  the  initial  critical  radius  of  the  system.  Each  bubble  later  grows 
or  shrinks  depending  upon  the  current  value  of  R^,,  which  depends  on  the  history  of 
bubble  size  changes  in  the  system.  Ultimately,  all  of  the  bubbles  collapse  or  else  only 
one  bubble  remains. 

Unless  specially  noted,  the  example  simulations  to  be  discussed  were  computed  at  the 
following  conditions: 

ambient  pressure  Q  =  1  ATA; 

temperature  T  =  293.16  °K  (20  °C); 

surface  tension  a  =  72.75  dyn/cm  (its  value  at  20  °C  for  the  air/H20  system); 

vapor  pressure  of  water  PH20  =  0.02307  atm  (its  value  at  20  °C); 

volume  of  the  condensed  phase  VL  =  0.015  cm3; 

diffusivity  =  3  •  10'5  cm/sec,  which  is  its  approximate  value  for  nitrogen  in  H20 
at  20°C7; 

partition  coefficient  K  =  6.34  •  10‘7  mol/(cm3-atm),  which  is  its  approximate  value 

for  nitrogen  in  H20  at  20  “C8. 
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The  system  volume  VL  also  will  be  seen  to  be  critical  to  this  analysis.  Its  value  of 
0.015  cm3  was  chosen  because  this  is  the  volume  of  a  membrane  of  dimensions  1  cm  x  1 
cm  x  150 Mm,  which  are  roughly  the  expected  dimensions  of  a  calibration  standard. 

Figure  3  shows  simulation  results  for  a  case  in  which  five  relatively  small  bubbles  are 
posited  in  a  condensed  phase  originally  containing  a  partial  pressure  of  1  atmosphere 
absolute  (1  AT  A)  of  dissolved  gas.  Rcrit  in  this  system  initially  is  60  Mm  (the  critical 
diameter  is  120  m  m).  All  of  the  bubbles  eventually  collapse,  the  smaller  ones 
disappearing  first,  the  partial  pressure  of  dissolved  gas  eventually  rises  to  1.004  ATA  but 
this  is  not  high  enough  to  prevent  the  largest  bubble  from  collapsing,  i.e.,  to  bring 
Rcrjt  <  R  for  the  largest  bubble. 

In  the  second  example,  shown  in  Figure  4,  all  conditions  are  the  same  as  in  Figure  3 
except  for  the  initial  sizes  of  the  bubbles.  Two  of  the  bubbles  initially  are  larger  than 
the  critical  diameter  of  120  Mm.  An  inspection  of  the  numerical  results  from  which  this 
plot  was  made  shows  that  the  bubble  that  initially  is  180  Mm  in  diameter  grows  slowly  for 
23  h  and  then  begins  to  shrink  very  slowly  as  F50  is  depleted  by  the  largest  bubble,  which 
will  continue  its  slow  monotonic  growth  as  it  asymptotically  approaches  equilibrium  with 
the  rest  of  the  closed  system. 

A  pseudo-stable  state  is  reached  in  Figure  4  because  two  of  the  bubbles  are  large 
enough  to  significantly  influence  the  amount  of  dissolved  gas.  Growth  of  a  sufficiently 
large  bubble  always  reduces  the  driving  force  for  diffusion  to  that  bubble  by  depleting 
P°°,  and  shrinkage  of  a  sufficiently  large  bubble  always  reduces  the  driving  force  for  that 
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Figure  3:  Bubble  size  changes;  "Small"  bubbles, 

initial  pcrtial  pressure  of  dissolved  gas  =  t  ATA 
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bubble  by  augmenting  P  00 .  Sufficiently  small  bubbles  won’t  affect  P  00  enough  to 
stabilize  their  sizes. 

In  Figure  5,  we  see  an  instance  of  a  substantial  initial  cverpressure  in  the  condensed 
phase.  The  initial  P 00  is  3  ATA,  giving  a  critical  diameter  of  1.4  pm.  Although  it  is 
difficult  to  see  on  the  plot,  the  bubbles  initially  are  very  different  in  size,  ranging  from 
5  pm  to  200  pm  diameter.  They  all  grow  rapidly  at  first;  the  smaller  ones  begin 
shrinking  as  P 00  declines.  At  time  =  5.8  h,  P 00  has  been  reduced  to  1  ATA  by  the 
bubble  growth.  Eventually,  only  one  bubble  will  remain.  But,  after  about  6  h  the 
bubbles  are  almost  static  in  size  over  many  hours.  This  is  because  the  driving  force  for 
mass  diffusion  is  small  since  P-P  00  and  the  surface  tension  term  2 a /F„  is  small  (see 
Equation  [15]). 

It  is  evident  from  studying  Figures  4  and  5  that  bubbles  having  diameters  much 
greater  than  100  Mm  remain  fairly  static  in  size  over  several  hours  when  the  overpressure 
is  small  and  the  surface  tension  is  that  of  the  air/H20  system.  That  is,  they  are  pseudo¬ 
stable  for  our  purposes.  Bubbles  much  smaller  than  100  pm  cannot  be  stabilized  at  a 
surface  tension  of  73  dyn/cm. 

Figure  6  illustrates  a  simulation  done  with  the  surface  tension  set  at  a  =  5  dyn/cm, 
which  may  be  possible  to  achieve  using  surfactants.  The  initial  value  of  P 00  is  1  ATA,  so 
that  a  slight  ‘overpressure’  exists  only  because  of  the  contribution  of  water  vapcr.  Still, 
the  critical  diameter  is  only  8.6  p m  because  of  the  low  surface  tension.  None  of  the 
bubbles  is  larger  than  60  pm  diameter  at  the  start.  The  results  are  qualitatively  similar 
to  those  in  Figure  5  in  that  we  eventually  obtain  a  pseudo-stable  system  after  an  initial 
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period  of  relatively  rapid  changes.  The  stability  again  occurs  because  of  the  small 
driving  force,  with  PP  00  and  2o/R  being  small.  In  this  hypothetical  system  we  have 
succeeded  ill  stabilizing  bubbles  smaller  than  200 

Can  bubbles  smaller  than  20  /urn  be  stabilized  if  both  the  surface  tension  and  the 


overpressure  are  extremely  small?  In  Figure  7,  P 00  is  0.99  ATA  initially,  the  surface 
tension  is  5  dyn/cm,  so  the  critical  diameter  is  15  *i  m.  The  bubbles  are  initially  5,  10,  15, 
20,  and  30  /on  in  diameter.  The  bubbles  that  begin  life  with  diameters  of  15  nm  or  less 
still  are  not  stabilized:  they  all  collapse  within  36  h.  The  larger  bubbles,  which  grow, 
will  not  stabilize  until  they  are  large  enough  to  deplete  the  dissolved  gas  surrounding 
them;  this  does  not  happen  within  the  first  100  h  of  the  simulation,  by  which  time  they 
both  are  larger  than  90  /im  in  diameter.  It  does  not  seem  possible  to  maintain  stable 
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although  it  certainly  would  be  possible  in  a  much  smaller  system. 
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Figure  7:  Bubble  size  changes;  "Small"  bubbles, 

Initial  partial  pressure  of  dissolved  gas  =  0.99  ATA, 
[  surface  tension  =  5  dyn/cm 
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VI.  NOISE  MEASUREMENTS 


We  attempted  to  quantify  the  amount  of  electronic  noise  generated  by  the  system’s 
various  noise  sources.  Of  primary  interest  is  the  amount  of  noise  present  when  the 
divide -by- 2  circuit  is  in  use,  because  this  is  the  output  from  which  we  hope  to  deduce 
bubble  sizes  by  analyzing  the  2nd  harmonic  component  of  the  sound  radiated  from 
bubbles.  Because  the  amplitude  of  harmonic  distortion  products  normally  depends  on 
the  amplitude  of  the  input  signal,  the  2nd  harmonic  noise  should  increase  with  the 
amplitude  of  the  input  signal.  This  is  borne  out  by  the  data. 

A.  Approach  and  methods 

In  the  following  discussion,  we  report  the  ratio,  in  decibels  (dB),  of  V  2h  to  Vmaxfund  for 
various  input  voltage  levels,  where  V  2h  is  the  RMS  (root  mean  square)  voltage 
amplitude  of  the  2nd  harmonic  distortion  and  Vm„fund  is  the  maximum  possible  RMS 
amplitude  voltage  of  the  fundamental  component  of  the  output  signal.  Vmaxfjnd  was 
measured  by  reflecting  the  transmitted  signal  off  of  the  surface  of  a  lead  brick  with  the 
divide-by-2  circuit  bypassed;  lead  is  an  excellent  sonic  reflector  because  of  the  mismatch 
in  mechanical  impedance  between  lead  and  water,  and  the  observed  output  signal 
amplitude  therefore  is  taken  to  be  maximized  in  this  experiment  for  a  given  input 
voltage  amplitude. 

V  2h  was  measured  by  reflecting  the  transmitted  signal  off  of  a  slab  of  Wall-Gone 
rubber  sound  absorber  with  the  divide-by-2  circuit  in  use  and  no  bubbles  apparent. 

Using  Wall-Gone  as  the  acoustical  target  is  more  realistic  as  a  physiological  analogue 
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than  measuring  noise  with  a  hard  target  in  place.  Although  ideally  we  might  expect  the 
signal  in  this  configuration  to  be  immeasurably  small,  in  fact  we  observe  a  substantial 
2nd  harmonic  component  generated  by  2nd  harmonic  distortion  in  the  front  end 
components  (consisting  of  the  transducers,  the  divide-by-2  circuit,  a  power  amplifier  in 
line  between  the  signal  generator  and  the  transmitting  transducer,  and  a  mixer),  and  by 
nonlinear  vibration  of  the  Wall-Gone  rubber.  On  the  other  hand,  noise  levels  are 
considerably  lower  with  Wall-Gone  as  the  target  than  with  a  hard  target,  because  less 
energy  is  reflected  to  the  receiving  transducer  from  the  Wall-Gone  than  from  a  hard 
target. 

The  scheme  for  characterizing  noise  levels  using  the  ratio  V  2h/Vmaxfund  was  chosen 
because  the  results  are  relevant  to  the  system’s  S/N  ratio  for  measurements  of  gas 
bubbles.  Although  the  reader  probably  finds  this  elevance  mysterious  now,  we  will 
explain  it  in  section  VIII,  ‘Estimating  the  Signal/Noise  Ratio’.  Note  for  now  that  V  2h, 
not  Vmaxfund,  is  considered  the  magnitude  of  noise. 

All  mea.  ements  were  made  under  the  following  conditions: 

'range"  mode  in  use 
sweep  time  =  0.1  sec 

span  =  4  kHz  (therefore,  time  record  length  =  0.1  sec) 
sweep  start  frequency  =  0.2  MHz 
sweep  st^'  frcqu  -  -  5  MHz 
windowing  =  flattop 

signal  averaging:  20  time  records  averaged 
target  distance  -  3.25  cm 
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Recall  from  the  Introduction  that,  in  the  range  mode,  the  output  signal  appears  as  a 
spectrum  of  voltage  amplitude  versus  offset  frequency  A  f,  where  A  f  is  proportional  to  the 
distance  from  the  transducer  head  to  the  target.  In  the  experiments  of  this  section  the 
output  spectrum  contained  a  voltage  spike  at  a  frequency  Af-2.1  kHz,  which  corresponds 
to  a  distance  of  3.25  cm  in  water.  The  reported  voltage  was  read  at  the  maximum  of  the 
peak  at  2.1  kHz.  Noise  at  other  frequencies  was  negligible  compared  with  the  noise  at 
2.1  kHz. 

B.  Results 

Table  1  shows  the  noise  level  measured  over  a  five-fold  range  of  input  voltage  levels. 
Inspection  of  the  rightmost  column  in  the  table  reveals  that  the  noise  level  is  essentially 
constant  at  — 37  dB  relative  to  the  maximum  output  signal  level  at  input  voltages  greater 
than  2  V. 

As  an  incidental  (but  important)  observation,  we  see  from  the  entries  in  the  third 
column  that  the  system’s  frequency  response  (i.e.,  output  voltage  at  the  fundamental 
frequency  -f  input  voltage)  is  nearly  linear  over  the  input  voltage  range  of  2.5  to  4  V 
peak-to-peak  from  the  signal  generator,  suggesting  that  this  is  an  advantageous  operating 
range.  From  additional  measurements  (not  shown)  it  has  been  observed  that  the  linear 
frequency  response  occurs  over  the  same  range  of  input  amplitudes  even  when  Wall- 
Gone  is  the  acoustical  target.  The  substitution  of  Wall-Gone  for  lead  attenuates  the 
output  signal  by  94%  consistently,  regardless  of  input  amplitude.  It  is  suggested 
therefore  that  the  appearance  of  nonlinear  frequency  response  is  a  much  stronger 
function  of  Vj(1  than  of  VQUt. 
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TABLE  1:  System  Noise 


voltage  setting  on 


3325A  signal  generator 

y  2h 

y  fund 

G 

> 

X 

> 

1  V  peak-to-peak 

47.42  m  V  RMS 

0.793  mV  RMS 

-31  dB 

1.5 

90.78 

1.523 

-35 

2 

136.1 

2.288 

-36 

2.5 

181.0 

3.032 

-36 

3 

229.6 

3.826 

-37 

3.5 

268.7 

4.43  v 

-38 

4 

301.6 

4.969 

-37 

4.5 

332.1 

5.398 

•37 

5 

357.3 

5.781 

-37 
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VII.  THEORY  OF  VIBRATING  BUBBLES 


A.  Approximate  analytic  solutions  via  the  perturbation  method 

The  motion  of  a  free,  spherical  gas  bubble  in  an  incompressible  condensed  phase  has 
been  described  by  solving  the  equation  of  motion,  i.e.,  the  momentum  balance,  for  the 
condensed  phase.  This  equation  is  found  to  be  an  ordinary  nonlinear  differential 
equation.  Since  the  system  is  spherically  symmetric  (given  a  homogeneous,  isotropic 
medium)  the  only  spatial  dimension  of  importance  is  the  radial  coordinate  r.  If  the 
equation  is  written  at  the  outer  bubble  surface,  we  obtain9,10 

RR"  +  (3/2)(R')2  =  1/P  J  Pin(t)  -  2o/R  -  tR  ]  (19) 

where  R  =  R(t.)  is  the  bubble  radius  at  time  t,  cm; 

R'  =  dR/dt; 

R'  -  d2R/dt2; 

pL  =  density,  g/cm3; 

Pjn(t)  =  the  pressure  at  the  inner  surface  of  the  bubble  at  time  t,  atm; 

o  =  surface  tension,  dyn/cm; 

t  R  =  radial  stress  at  the  outer  bubble  surface,  dyn/cm2. 

The  term  tR  arises  when  we  apply  the  boundary  condition  that  radial  stresses  are 
continuous  across  the  bubble  surface.  In  section  C  we  show  how  this  stress  term  is 
rewritten  in  terms  of  strain  or  rate  of  strain,  depending  on  what  constitutive  equation  is 
valid  for  the  condensed  phase.  For  now,  let  us  assume  that  the  bubble  resides  in  a 


Newtonian  fluid.  In  section  C  we  will  show  that  this  leads  us  to  calculate  xR  =  Pout  + 
4fiR'/R,  where  Pou‘  is  the  hydrostatic  pressure  at  the  outer  surface  of  the  bubble. 

Prosperetti 10,11  found  an  approximate,  analytic  solution  to  Equation  [19]  for  a 
sinusoidally  oscillating  ambient  pressure,  that  is  for  a  bubble  in  an  externally  applied 
sound  field.  The  mathematics  are  quite  complicated  and  we  will  only  summarize  some 
of  the  major  features  here.  For  more  details  see  the  Appendix  of  reference  1.  For  a 
reasonably  complete  discussion  of  the  mathematics  one  must  look  to  Prosperetti’s 
papers10,11. 

Prosperetti  replaces  Pout(t)  with  P0(l-ecos(Q  t)),  i.e.  P0  is  the  mean  ambient  pressure, 
e  is  the  amplitude  of  sound  waves  imposed  on  the  bubble,  and  the  frequency  of  that 
sound  is  Q .  The  internal  pressure  Pin(t)  is  represented  by  a  polytropic  expression 
Pjn  =  Pjneq  (R0/R(t))3‘  -  4*ithR'/Ro>  where  Ro  is  the  equilibrium  radius,  k  is  the 
"effective  polytropic  exponent",  Pln  is  the  internal  pressure  corresponding  to  R=R<t,  and 
uth  is  a  so-called  ‘thermal  viscosity’  introduced  as  an  ad  hoc  corrective  term  to 
compensate  for  the  polytropic  expression’s  inability  to  account  for  phase  shifts  between 
pressure  and  temperature.  Prosperetti  also  defines  an  ‘acoustic  viscosity’  /xac  that  arises 
from  momentum  transfer  from  the  bubble  to  the  liquid.  The  thermal  and  acoustic 
viscosities  have  the  same  units  as  a  "real"  viscosity  and  each  is  a  "damping"  term, 
meaning  that  it  appears  in  the  momentum  balance  as  a  multiplier  to  the  velocity  term 
R'.  Accordingly,  Prosperetti  combines  all  three  viscosities  into  an  "effective  viscosity" 
Men-  =  +  Mth  +  Mac-  Equation  [19]  now  becomes 
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RR’  +  (3/2)(R  )2  = 

1/pJ  Pjn,eq  (V^O)"  -  P0(l-€COS(Qt))  -  2o /R  -  4/,cffR '/R  ]  [20] 

Expressions  for  the  polytropic  exponent,  thermal  viscosity,  and  acoustic  viscos  ity  all  are 
calculated  from  a  separate  analytic  solution  to  the  system  of  six  partial  differential 
equations  that  comprise  the  linearized  equations  of  continuity,  motion,  and  energy  in 
both  phases.  The  condensed  phase  is  considered  compressible  in  the  linearized 
governing  equations,  even  though  the  nonlinear  Equations  [19]  and  [20]  are  written  for 
incompressible  liquids.  (In  fact,  the  acoustic  damping  term  is  zero  for  an  incompressible 
medium.)  linearized  equations  are  valid  for  small  perturbations  from  equilibrium  --  in 
this  case,  for  small  amplitude  vibrations  --  but  keep  in  mind  that  no  information  on 
harmonics,  subharmonics,  or  phase  shifts  is  contained  in  a  linearized  treatment  of  a 
vibrating  system.  The  mechanics  of  solving  this  system  of  equations  are  exceptionally 
complicated. 

Next,  the  substitution  R(t)  =  R0  (1  +  x(t))  is  made  into  equation  [20],  The 
dimensionless  variable  Y  is  seen  to  be  a  fractional  perturbation  of  the  bubble  radius 
R(t)  from  its  equilibrium  value  R^  The  problem  is  made  tractable  by  discarding  ail  x,  x’, 
and  x"  terms  that  are  raised  to  powers  of  4  or  greater,  but  note  that  retaining  the  2nd 
and  3rd  power  terms  retains  the  essential  effect  of  the  nonlinearity  (provided  x  is  small). 

Finally,  algebraic  solutions  to  the  resulting  equation  are  obtained  by  a  "perturbation" 
technique  first  advanced  by  Krylov,  Bogolyubov,  and  Mitropolsky12.  Each  solution  is 
valid  when  the  forcing  frequency  Q  is  ‘near’  one  of  the  system’s  resonance  frequencies. 
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We  arc  particularly  interested  in  the  2nd  harmonic  component  of  the  solutions,  which  is 
found  by  retaining  only  the  terms  that  arc  of  frequency  2Q,  since  that  component  will  be 
measured  by  the  bubble  detector.  It  turns  out  that  the  mathematical  system  (and 
presumably  the  physical  system  likewise,  if  its  mathematical  description  is  accurate)  is 
‘catastrophic’:  as  the  forcing  frequency  is  swept  through  the  various  resonance 
frequencies  there  are  discontinuous  jumps  between  stable  solutions,  and  there  is 
hysteresis  behavior  that  causes  the  location  of  these  jumps  to  depend  on  whether  the 
forcing  frequency  is  being  swept  upward  or  downward. 

The  procedure  summarized  above  yields  an  algebraic  expression  for  the  amplitude  of 
the  oscillations  of  R  as  a  function  of  forcing  frequency  and  all  other  relevant  parameters. 

This  is  readily  converted  to  the  amplitude  of  the  radiated  pressure  wave  corresponding 
to  it,  denoted  l|PnKll|,  by  solving  the  linearized  *wave  equation’  in  the  liquid  (in  fluid 
mechanics  a  wave  equation  is  simply  the  equation  of  motion  written  for  a  compressible 
medium  in  which  a  wave  is  being  propagated).  If  desired,  only  the  2nd  harmonic 
component  of  the  radial  oscillations  can  be  included,  and  we  designate  the  pressure 
amplitude  of  that  component  as  ||Prad2h||. 

The  fact  that  the  values  of  the  parameters  ‘k’  and  7a’  (used  below)  are  estimated 
using  first  order  expansions  (linearizations)  of  the  terms  in  the  mass,  energy,  and 
momentum  balances  means  that  they  are  valid  only  for  a  linear  system.  Their  usage 
when  analyzing  this  nonlinear  system  is  an  ad  hoc  approach  necessitated  by  the  difficulty 
of  the  mathematics. 

Clearly,  the  analytic  solution  incorporates  several  questionable  assumptions  and 
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computational  simplifications.  Its  value  is  that  it  can  be  computed  rapidly  enough  to 
make  it  useful  in  the  iterative  data-fitting  routine  mentioned  in  the  Introduction  to  this 
report,  with  which  the  equilibrium  bubble  radius  will  be  estimated  for  an  unknown 
bubble  by  comparing  observed  resonance  frequencies  with  the  predicted  values. 

B.  Numerical  solutions  via  the  Galerkin  Spectral  Method 

Given  our  reservations  about  the  reliability  of  the  analytic  solution  discussed  in 
section  A,  we  would  like  to  have  a  means  of  numerically  solving  the  equations  governing 
a  vibrating  bubble  without  resorting  to  so  many  simplifying  assumptions.  This  new 
solution  could  be  used  to  evaluate  the  accuracy  cf  the  analytic  solution  in  the  absence  of 
direct  experimental  verification,  which  is  unlikely  to  appear  any  time  soon.  We  are 
willing  to  trade  computational  economy  for  accuracy  with  the  numerical  solution. 

Kamath  and  Prosperetti13  used  a  Galerkin  spectral  method  to  solve  the  governing 
equations  and  found  significant  deviations  of  the  analytic  solution  from  the  numerical 
results  when  the  driving  amplitude  was  high  enough.  As  will  be  discussed  shortly,  the 
mathematical  description  that  Prosperetti  evaluated  using  the  Galerkin  method  is  in 
some  ways  more  sophisticated,  in  other  ways  more  simplified  than  the  description  he 
evaluated  earlier  using  a  perturbation  technique.  This  confuses  a  comparison  of  the 
results  from  the  two  solutions.  The  reader  should  look  to  reference  14  for  an  overview 
of  Galerkin  techniques.  Here,  we  will  review  the  particulars  of  how  Prosperetti  (and  we) 
used  the  Galerkin  spectral  method  on  our  system. 

The  Galerkin  method  is  used  to  reduce  the  partial  differential  equations  to  a  system 
of  ordinary  differential  equations,  which  subsequently  are  solved  using  the  Gear  method. 
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The  result  consists  of  the  value  of  R,  the  bubble  radius,  as  a  function  of  time  at  some 
forcing  frequency  Q .  This  time  domain  information  can  be  converted  to  an  amplitude 
spectrum  via  Fourier  transformation.  To  get  the  steady  state  solution  one  integrates 
until  steady  state  oscillations  are  attained,  then  Fourier  transforms  the  time  domain  data 
from  only  the  last  oscillation  cycle.  If  one  is  interested  in  the  amplitude  of  only  the  2nd 
harmonic  component  as  a  function  of  Q ,  then  the  amplitude  at  frequency  2Q  is  extracted 
from  the  spectrum.  Once  the  steady  state  radial  amplitude  at  frequency  2Q  has  been 
computed  over  a  range  of  Q ’s,  we  can  construct  the  spectrum  of  amplitude  versus  Q . 
Also,  the  amplitude  of  R  at  frequency  2Q  can  be  converted  to  ||  Prad2h||  by  solving  the 
aforementioned  linearized  wave  equation  in  the  condensed  phase.  In  sections  D  and  E, 
we  will  present  results  in  an  alternative  form  that  does  not  involve  Fourier 
transformation:  we  show  the  maximum  bubble  radius  attained  during  a  steady  state 
oscillation. 

In  Prosperetli,  Crum,  and  Commander,15  simplified  governing  equations  are  derived 
from  the  conservation  equations  for  mass,  momentum,  and  energy  for  both  the  gas  in  the 
bubble  and  the  liquid  outside  the  bubble,  along  with  boundary  conditions  at  the  gas- 
liquid  inteiface  (i.e.,  the  bubble  wall).  Based  on  several  reasonable  assumptions,  and 
several  order  of  magnitude  comparisons  that  show  certain  quantities  to  be  negligible,  the 
authors  are  able  to  reduce  the  original  system  of  six  partial  differential  equations  (i.e., 
the  conservation  equations)  to  one  partial  differential  equation  and  two  ordinaiy 
differemial  equations. 

Two  basic  assumptions  are  that  the  gas  is  perfect  and  the  bubble  maintains  a 
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spherical  shape.  An  order  of  magnitude  argument  shows  that,  near  the  main  resonance 
frequency  of  a  bubble  in  forced  oscillation,  the  Mach  number  MB  of  the  bubble  wall, 

referred  to  the  speed  of  sound  in  the  gas,  is  extremely  small;  i.e.  MB  -  R/cc  is  close  to 
zero.  Using  this  approximation  it  can  be  shown13  that  the  radial  gradient  in  pressure 
inside  the  bubble  is  negligible  compared  with  the  temporal  variation  of  pressure;  i.e.,  the 
pressure  is  almost  uniform  spatially  within  the  bubble.  An  estimate  of  viscous  shear  in 
the  gas  shows  it  to  be  negligible.  Based  on  these  results  it  can  be  shown  that  the 
momentum  balance  equation  for  the  gas  reduces  to  the  simple  statement  that  gas 
pressure  is  a  function  only  of  time. 

Some  further  assumptions,  which  involve  the  bubble  wall  temperature  and  vapor 
effects,  are  less  important,  because  the  system  equations  would  still  be  tractable  without 
them.  An  order  of  magnitude  estimate  shows  that,  owing  to  the  greater  thermal  inertia 
of  liquid  than  gas,  the  temperature  at  the  bubble  surface  is  almost  unperturbed  from  the 
ambient  temperature  provided  there  is  not  too  much  heat  being  alternately  consumed 
and  released  by  the  evaporation  and  condensation  of  vapor.  The  authors  ensure  the 
validity  of  this  approximation  in  their  mathematical  construct  by  neglecting  vapor 
altogether.  In  the  physical  world,  they  suggest  (based  on  their  order  of  magnitude 
estimates)  that  the  approximation  is  valid  for  bubbles  in  H20  at  temperatures  up  to 
about  50  °C.  In  other  words,  the  liquid  must  be  sufficiently  cool.  With  the  use  of  this 
approximation  the  equation  of  energy  in  the  liquid  phase  can  be  dispensed  with; 
otherwise,  it  would  be  needed  for  determining  the  boundary  condition  on  the  energy 
equation  in  the  gas  phase. 
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In  addition,  both  slow  and  fast  diffusion  of  the  gas  in  and  out  of  the  bubble  are 
known  to  have  negligible  dynamical  effects.  Thus,  diffusion  is  ignored  entirely  and  we 
assume  the  bubble  boundary  is  impervious  to  the  gas.  Fast  diffusion  refers  to  diffusion 
with  time  scale  the  order  of  one  period  of  the  oscillation;  slow  diffusion  (often  called 
rectified  diffusion)  refers  to  diffusion  with  a  time  scale  of,  for  example,  thousands  of 
periods  in  which  the  total  gas  contents  and/or  the  equilibrium  radius  of  the  bubble  may 
change. 

Under  the  stated  assumptions,  the  system  of  equations  to  be  solved  becomes: 
a)  Radial  dynamics  equation: 

We  use  Keller’s  equation  for  a  gas  bubble  in  a  compressible  liquid. 

(1  -R/c)RR  +  3/2(1  -R/  (3c))  R2  - 

l/pL(l  *  R/c  *  (R/c)  djdt)p3  -  1/Pi(l+  R/c)pA(t  -  R/c)  [21] 

where  R  -  R(t)  is  the  bubble  radius  at  time  t,  p  L  is  the  density  of  the  liquid,  ’c'  is  the 
speed  of  sound  in  the  liquid,  nB  is  the  pressure  in  the  liquid  just  outside  the  bubble,  and 
pA(t)  is  the  ambient  pressure  (in  radial  equation  PA  is  evaluated  at  time  =  t  +  R/c). 
Keller’s  equation  is  valid  when  the  Mach  number  of  the  bubble  wall  is  sufficiently  small 
so  that  quantities  proportional  to  its  square  may  be  neglected.  We  assume 
pA(t)  m  p„  (1  -  e -suifur)) .  A  balance  of  the  normal  stresses  across  the  bubble  surface 

yields: 
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Ps  “  PgW  "  P  -  2a/R  -  4  Vl(r/r) 


[22] 


where  -  liquid  viscosity, 
b)  Pressure  equation: 

This  equation  relates  the  internal  pressure  to  the  temperature  gradient.  Since  the 
internal  pressure  p  is  assumed  to  be  uniform  throughout  the  bubble;  p  is  independent  of 
r;  it  depends  only  on ’t’. 

p-3/tf(  (y  -  l)K‘(dTfdr)\x  -  ypR)  [23] 


c)  Temperature  equation: 

This  is  also  called  the  energy  equation. 

YP/((Y  "  1)7)  [dTfdt  +  1/(y  p)  ((Y  -  1)  K  dT/dr  -  U3rp)d7]dr]  -  p  -  V  •  (K  VT)  | 

[24] 

where  K(T)  is  the  thermal  conductivity  of  air.  Based  on  data  for 

200  °K  <  T  <  3000  °K,  Kamath  and  Prosperetti  suggest  using  the  linear  function 

K(T )  -  5.528  T  +  1165. 

Boundary  conditions  for  system: 

The  temperature  on  the  bubble  surface  equals  the  constant  exterior  temperature  Tm . 
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Initial  conditions  for  system: 

Bubble  is  at  rest  with  equilibrium  radius  R0 . 

Internal  pressure  has  value  p0  m  (2  o  /  R0 )  +  pm 

Temperature  is  constant  throughout  bubble;  its  value  is  Tm . 

Let  us  define  a  scaled  radial  variable  y  ■  r/R(t)  to  be  used  instead  of  r. 

A  dependent  variable  x  (y,t)  is  defined  by  x  ■  fT  K(G)  dO  where  K(T)  is  defined 

Jr. 

above. 

The  Galerkin  method  is  used  to  replace  spatial  partial  derivatives  (i.e.,  those 
involving  the  variable  *y’)  with  finite  linear  combinations  of  Chebyshev  polynomials 
Tk(y ) .  Specifically,  we  represent  x(y,t)  as  an  infinite  sum  of  shifted  Chebyshev 
polynomials  *t(y)  «  Tlk(y)  -  1 .  Let  the  coefficients  be  denoted  ak(t) .  Denote  by 

N 

x  N  the  series  truncated  after  N  terms;  xN  -  Y*  <**(0 

*-i 

The  pressure  and  temperature  equations  involve  spatial  derivatives.  Using  the 
expression  for  we  can  convert  these  spatial  derivatives  into  finite  linear  combinations 

N 

^2  c*(0  where  the  ek{t)  are  functions  of  the  ak(t) . 

k-  1 

With  these  techniques,  the  temperature  equation  is  reduced  to  a  system  of  ordinary 
differential  equations  which  can  be  integrated  in  time.  Following  Kamath  and 
Prosperetti’s  approach  we  used  a  four-term  Chebyshev  expansion.  The  system  must  be 
integrated  for  many  (usually  more  than  10)  cycles  before  the  system  reaches  steady-state 
oscillations.  For  certain  parameter  values  there  may  be  more  than  one  stable  solution;  if 
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so,  the  solution  that  arises  in  any  giv  :n  problem  will  depend  on  the  initial  conditions 
used  for  that  problem.  In  other  words,  the  bistability/hysteresis  effect,  well  established 
for  the  analytic  solution  with  resonance  frequencies  (see  the  preceding  section),  appears 
as  well  in  the  numerical  solution. 

C.  Derivation  of  the  normal  stress  at  the  surface  of  a  spherical  cavity 

The  boundary  condition  on  the  equation  of  motion  in  the  medium  surrounding  a 
bubble  (Equation  [19])  is  different  for  a  liquid  than  for  a  solid  because  different 
constitutive  equations  govern  the  mechanical  behavior  of  the  two  types  of  media. 

The  continuity  equation  for  an  incompressible  substance  with  spherical  symmetry  can 
be  written  either 


n  /iA  ar.-2™  \/Ar  =  n 

W  *  /  ~  V*  ’r// 

Is) 

JL 

or 

(l/r2)  3<rSir)/ar  =  0 

[25b] 

where  vr  =  radial  velocity  and  ur  =  radial  displacement.  Velocity  is  simply  the  time 
derivative  of  displacement:  vr  =  dur/c?t.  Integration  yields 


Vr  -  A,(t)/r 

[26a] 

or 

ur  =  A  2(t)/r 

[26b] 

where  A  j(t)  and  A  2(t)  are  integration  factors  to  be  evaluated  later.  The  boundary 
condition  on  the  continuity  equation  is  continuity  of  either  velocity  or  displacement  at 
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the  bi,  jble  surface: 


I 


Vrlr-R  =  dR/dt  =  R' 


urlr=R  ~  R  "  Re 


where  RE  =  bubble  radius  at  which  the  medium  surrounding  the  bubble  is  undcformed, 
which  in  general  is  neJ  the  same  as  the  equilibrium  radius  Ro.  It  follows  from  Equations 
[26]  and  [27]  that 


vr  =  R^'/r2 
ur  =  R2(R-RE)/r 


The  boundary  condition  on  the  equation  of  motion  is  continuity  of  the  radial  stresses 


across  the  phase  interface*: 


-Pjn  +  2o/R  =  tR 


where  x  R  is  the  racial  stress  at  r  =  R. 


We  arc  obeying  the  following  incomprehensible  arbitrary  convention  on  signs  of  stress  components, 
taken  almost  verbatim  from  reference  9: 

The  stress  r(j,  due  to  action  by  material  on  the  positive  side  of  the  svrface  on  the  material  on  the  negative 
side,  is  positive  if  the  line  of  action  is  along  positive  spatial  coordinate  x , . 

Conversely,  the  stress  exerted  from  the  negative  side  of  the  surface  on  the  material  on  the  positive  side  is 
positive  if  the  line  of  action  is  along  negative  x , . 
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Now  we  will  apply  the  constitutive  equations.  For  a  Newtonian  liquid,  the  stress 
teasor  components  are  linearly  related  to  Ihe  components  of  the  rate  of  strain  tensor: 


=  -P6- 


[30] 


where  P  =  hydrostatic  pressure,  atm; 

6  ij  =  Kronecker  delta,  defined  thus: 

=  1  when  i=j,  and  =  0  when  i*j; 

H  =  viscosity,  dyn/cm; 

A  =  component  of  the  rate  of  strain  tensor. 


The  components  of  the  linearized  rate  of  strain  tensor  are  linear  functions  of  velocity 
vector  gradients.  For  satisfying  the  boundary  condition  with  spherical  symmetry  we  need 
only  know  the  radial  normal  stress  x  n.  It  is  found  from  the  radial  normal  component  of 


♦h-a  r- 
uiv  * 


^"tv  of  strain,  WaIic!*  m  Sjplicrica.1  cooxuIuai^d  is  ^iv&n  uy 


Arr  =  2  d\Jdr 


[31] 


And  so  the  boundary  condition  equation  [29]  is  solved  by  xR  =  tn|rsR,  which  is  readily 
evaluated;  the  boundary  condition  for  a  Newtonian  liquid  becomes 

P0Ut  +  2o/R  =  Pta  -  4#iR'/R  [32] 

For  a  linearly  elastic  solid,  by  definition  the  stress  tensor  components  are  linearly  related 
to  the  components  of  the  strain  tensor16: 


Ty  —  A.  («]] +632+633)6  y  +  2^eij  [^3] 

where  k  and  n  are  the  Lame  constants; 
e  =  small  strain  tensor  component; 

6 11^22'  ar*d  c33  are  the  normal  strain  components. 

(By  longstanding  convention,  the  second  Lame  constant  unfortunately  is  denoted  by  the 
same  Greek  letter  as  is  the  viscosity  of  a  liquid.)  Equation  [33]  does  not  appear  at  first 
to  contain  the  hydrostatic  pressure,  but  we  can  make  the  substitution16 

P  =  -  B  (en +6^+633) 

=  -(X  +2  )  (£11+622+633)  [34] 

where  B  =  bulk  modulus,  dyn/cm2. 

Equation  [33]  now  becomes 

ty  =  [*P-(3p/3X611+e22+e33)]6ij  +  2^  [35] 

The  continuity  equation  tells  us  that  in  an  incompressible  system  the  sum  of  the  normal 
strain  components  (6H+622+633)  is  zero,  so  we  can  simplify  the  constitutive  relation  to 

ry  =  -P&y  +  2^6y  [36] 
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In  linearized  elastic  theory  the  strain  tensor  components  are  linear  functions  of 
displacement  vector  gradients.  For  this  spherically  symmetric  system  we  need  to  obtain 
the  radial  normal  stress  from  the  radial  normal  strain,  which  in  spherical  coordinates  is 
given  by 

err  =  dur/a  r  [37] 

And  so  the  boundary  condition  again  is  solved  by  tr  =  xrr|r=R.  This  radial  normal 
stress  component  is  given  by 


^rrir-R  =  P  -  ^  (R-R^/R  [38] 

at  the  bubble  surface.  Finally,  let  us  replace  the  Lame  constant  with  the  more  familiar 
Young’s  modulus  E.  For  an  incompressible  solid  p  =  E/3.  Our  final  expression  for  the 
linearized  elastic  boundary  condition  is  therefore 

Pout  +  2o/R  =  Pin  -  (4E/3)(R-Re)/R  [39] 

Linearized  elastic  theory  is  valid  for  small  deformations.  At  larger  strains  we  require 
the  nonlinear  theory,  in  which  expressions  for  the  strain  components  are  calculated  from 
gradients  in  the  displacement  vector  without  discarding  terms  of  order  2  and  higher.  The 
mathematics  involved  in  computing  our  boundary  condition  for  this  kind  of  system  are 
well  beyond  the  scope  of  this  report.  Fortunately,  the  result  has  been  published17: 
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tJ,.R  =  -P-(1/6)E  [5  -  4(RU/R)  -  (R^R)4] 


[40] 


The  boundary  condition  on  the  equation  of  motion  for  large  elastic  deformations  is  given 
by 

Pout  +  2o/R  =  Pin  -  (1/6)E  [5  -  4(Re/R)  -  (RE/R)4  [41] 

This  boundary  condition  approaches  the  one  for  a  linear  elastic  solid  (equation  [39])  in 
the  limit  as  R  -  RE. 

So,  the  appropriate  boundary  condition  is  substituted  into  the  equation  of  motion 
(Equation  [19])  and  that  equation  is  solved  either  analytically  or  numerically  to  yield  a 
description  of  an  oscillating  gas  bubble.  To  obtain  the  analytic  solution  it  is  necessary  to 
assume  linear  elasticity  when  modelling  oscillations  in  a  solid,  so  Equation  [39]  must  be 
used  to  calculate  the  elastic  stress  on  the  bubble  wall  arising  from  oscillations  of  the 
radius  about  its  static  value  Rq.  However,  the  nonlinear  Equation  [40]  is  used  to 
compute  the  elastic  stress  imposed  on  the  wall  of  the  static  bubble,  i.e.  the  stress  arising 
from  the  difference  between  Rq  and  RE.  Thus,  in  the  analytic  solution  we  have  a  linear 
oscillatory  stress  superimposed  on  a  nonlinear  static  stress.  When  using  the  Galerkin 
method  the  nonlinear  Equation  [40]  is  used  for  all  calculations  of  elastic  stress. 

There  is  one  potential  complication  when  figuring  the  boundary  condition  for  a 
vibrating  bubble  in  an  elastic  solid.  The  following  relations  are  true  for  the  strain  and 
elastic  stress  in  the  solid  surrounding  the  bubble: 


58 


and 


for  R  >  £„  <  0  and  <  0, 

for  Rs  Re  6^  =  0  and  in  =  0, 


[42a] 

[42b] 


where  Re  --  bubble  radius  at  which  the  medium  surrounding  the  bubble  is  undeformcd. 
(The  strain  and  stress  are  negative  when  R  >  RE  because  of  our  arbitrary  sign 
convention.)  Equation  [42]  says  that  the  bubble  pushes  against  the  solid  when  R  >  RE, 
but  doesn’t  pull  on  it  when  R  <  RE;  that  is,  the  solid  Goes  not  adhere  to  the  bubble  wall. 
Consequently,  the  boundary  condition  is  discontinuous  at  R  =  RE.  If  a  bubble  is  made 
to  oscillate  in  size  such  that  its  radius  is  greater  than  RE  during  expansions  and  less  than 
Re  during  contractions,  then  its  boundary  condition  is  discontinuous  in  time.  Describing 
such  a  system  mathematically  would  present  formidable  challenges.  We  therefore  limit 
ourselves  to  analyzing  systems  for  which  R  is  never  less  than  RE  during  vibration. 
Experimentally,  if  a  bubble  is  formed  by  cavitation  in  a  solid  then  RE  =  0  and  the 
discontinuous  boundary  condition  can  never  exist.  If,  on  the  other  hand,  the  solid  phase 
is  formed  by  solidifying  a  liquid  around  a  gas  bubble  that  is  at  its  equilibrium  radius 
then  we  expect  that  Re  =  Rq,  and  if  the  bubble  radius  is  made  to  oscillate  about  R0  (in 
response  to  sound  waves,  for  instance)  then  the  boundary  condition  is  discontinuous  in 
time. 

D.  Results:  Bubbles  in  a  viscous  liquid 

Using  FORTRAN  programs  that  implement  the  models  discussed  above,  we  have 
generated  the  graphs  shown  below.  In  all  these  graphs,  we  show  the  maximum  bubble 
radius  attained  during  an  oscillation  after  the  system  is  in  a  steady  state  mode.  The 
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maximum  value  is  given  as  a  function  of  the  driving  frequency;  the  latter  is  cxpiessed  as 
a  ratio  of  driving  frequency  to  natural  frequency  of  the  bubble. 

The  expression  ‘Galerkin  Method’  in  the  graph  headings  refers  to  the  program  that 
implements  the  Keller  equation  (see  above)  for  bubble  oscillations.  In  this  program  the 
differential  equation  modeling  heat  flow  across  the  bubble  wall  is  solved  numerically 
using  the  Galerkin  spectral  method.  The  expression  ‘Penurbation  Method’  refers  to  the 
program  which  implements  Rayleigh’s  equation  (with  the  addition  of  thermal  and 
acoustic  regarding  damping  terms  as  discussed  above,  and  a  natural  frequency  <a0  for  the 
bubble  which  depends  on  the  driving  frequency  <*>).  The  graphs  are  based  on 
approximate  analytic  solutions  proposed  by  Prosperetti10.  These  two  methods  are  used 
for  the  case  of  a  bubble  in  water  driven  by  a  sound  wave  with  a  small  amplitude  (0.3 
atm).  For  the  associated  graphs,  see  Figures  8  and  9.  Figures  10  and  11  are  graphs  for 
the  case  of  a  bubble  in  water  driven  by  a  sound  wave  with  a  larger  amplitude  (0.6  atm). 
For  this  larger  amplitude,  which  we  use  as  well  for  the  last  four  figures,  we  use  another 
program  based  on  Rayleigh’s  equation  without  the  enhancements  described  above.  The 
program,  which  we  refer  to  as  the  ‘Simple  Model’,  solves  the  differential  equation  using 
a  very  accurate  numerical  method.  Use  of  a  numerical  solution  makes  unnecessary 
certain  analytic  approximations  required  in  the  perturbation  approach.  The  greater 
numerical  accuracy  achieved  allows  the  simple  model  to  handle  large  driving  amplitudes 
which  lead  to  numerical  instabilities  using  the  perturbation  theory  approach. 
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Figure  8:  Maximum  Bubble  Radius  (re!,  to  Ro);  Galerkin  Method 
Bubble  in  H2O;  Ro  =  10^m  Driving  Amplitude  =  .3  atm. 


jure  9:  Maximum  Bubble  Radius  (rel.  to  Ro);  Perturbation  Method 
Bubble  in  H2O;  Ro  =  1 0^m  Driving  Amplitude  =  .3  atm. 


Figure  11:  Maximum  Bubble  Radius  (rel.  to  Ro);  Simple  Model 
Bubble  in  H20;  Ro  =  lO^m  Driving  Amplitude  =  .6  atm. 


The  last  four  figures  apply  to  the  case  of  an  elastic  solid  with  a  shear  modulus  of 
8  •  105  dyne/cm2.  This  is  a  low  value  for  the  shear  modulus  and  indicates  a  soft  elastic 
solid.  For  example,  this  shear  modulus  would  be  typical  of  a  fairly  high-concentration 
polyacrylamide  gel.18  Figure  12  (Galerkin  Method)  and  Figure  13  (Simple  Model) 
describe  the  case  of  a  bubble  for  which  Rq  =  Rc  (recall  that  the  latter  is  the  radius  at 
which  there  is  no  strain  in  the  solid).  This  means  that,  at  mechanical  equilibrium,  the 
solid  is  unstrained  and  the  bubble  wall  is  not  subject  to  elastic  stress.  We  use  a  common 
value  of  10  because  bubbles  of  this  size  are  of  interest.  For  Figure  14  (Galerkin 
Method)  and  Figure  15  (Simple  Model),  we  consider  the  case  of  a  bubble  for 
which  Ro  >  Re.  We  allow  for  large  deformations;  the  nonlinear  stress-strain  relation  is 
described  by  Equation  40.  Here,  R^,  =  10 /am,  RE  =  5  fim. 

In  general,  there  are  two  stable  solutions  to  both  Rayleigh’s  and  Keller’s  equation  for 
values  of  the  frequency  ratio  that  lie  in  the  resonance  bands.  The  perturbation  program 
will  always  find  both  these  solutions  since  it  has  analytic  expressions  for  them.  The 
situation  for  the  Galerkin  program  and  the  simple  model  program  is  more  involved.  The 
stable  solution  that  is  numerically  calculated  in  these  programs  depends  on  the  initial 
conditions;  in  general,  starting  with  the  bubble  at  rest  and  with  radius  Rj,  leads  to  the 
solution  with  the  smaller  amplitude.  Finding  initial  conditions  that  lead  to  another 
stable  solution  involves  some  numerical  experimentation.  For  these  differential 
equations  only  at  most  two  stable  solutions  have  been  found  to  exist  in  any  resonance 
band;  however  it  is  possible  others  exist.  For  the  case  described  in  Figure  13  we  were 


able  to  find  such  initial  conditions  (initial  radius  =  2  •  Rg  cm,  and  initial  velocity  =  Rg 
cm/sec);  therefc  e  we  display  both  solutions.  As  the  system  is  stiffened,  either  by 
increasing  the  shear  modulus  or  by  decreasing  R&  the  two  solutions  may  coalesce  into  a 
single  one.  One  has  to  consider  this  possibility  when  one  is  unable  to  find  more  than 
one  stable  solution. 
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Figure  12:  Maximum  Bubble  Radius  (rel.  to  Ro);  Galerkin  Method,  shear  =  8e5, 
Bubble  in  Elastic  Gel;  Ro  =  Re  =  10/um  Driving  Amplitude  =  .6  atm. 


Figure  13:  Moximum  Bubble  Radius  (rel.  to  Ro)i  Simple  Model 
Bubble  in  Elastic  Gel;  Ro  =  Re  =  10/M"  Driving  Amplitude  =  .6  atm 


Figure  14:  Maximum  Bubble  Radius  (rel.  to  Ro);  Galerkiri  Method,  shear  =  Se5, 
Bubble  in  Elastic  Gel;  Ro  =  lO^um  Re  =  5/xm  Driving  Amplitude  ~  ,6  atm. 


Figure  15:  Maximum  Bubble  Radius  (rel.  to  Ro);  Simple  Model,  shear  =  8e5, 
Bubble  in  Elastic  Gel;  Ro  =  10/xm  Re  =  5/xm  Driving  Amplitude  =  .6  atm. 


VIII.  ESTIMATING  THE  SIGNAL/NOISE  RATIO 


We  can  use  the  models  presented  in  section  VII  to  estimate  the  ratio  of  (the  sound 
pressure  amplitude  at  the  2nd  harmonic  frequency  radiated  by  a  spherical  bubble  to  the 
receiving  transducer]  to  (the  sound  pressure  amplitude  at  the  fundamental  frequency 
radiated  by  the  transmitting  transducer].  This  tells  us  the  "line  loss"  attendant  with 
measuring  a  bubble.  For  instance,  if  Pout  is  one-twentieth  of  Pin  then  the  line  loss  is 
95%  or  13  dB.  If  we  assume  that  the  bubble  detector’s  overall  frequency  response  is 
linear  (that  is,  Vout  is  a  linear  function  of  Vin)  then  for  this  example  Vout  would  be  13  dB 
below  its  maximum  attainable  value. 

In  section  VI,  ‘Noise  Measurements’,  we  reported  a  relative  noise  level  equal  to  the 
ratio  of  [the  2nd  harmonic  noise  level]  to  [the  maximum  attainable  level  of  the 
fundamental  of  the  output  signal].  Therefore,  by  inspection,  we  can  estimate  the  S/N 
ratio  of  the  system  by  dividing  the  line  loss  associated  with  a  given  bubble  by  the  relative 
noise  level  as  reported  in  section  VI.  In  this  section  we  discuss  the  estimated  S/N  ratio 
for  single  spherical  bubbles  of  various  sizes  in  water. 

The  line  loss  estimates  were  made  using  the  analytic  solution  to  the  governing 
equations.  The  line  loss  is  predicted  to  be  only  a  weak  function  of  the  amplitude  of  the 
driving  signal  within  the  range  of  0  to  0.5  atm  sound  pressure  amplitude.  The  analytic 
solution  is  valid  only  at  small  amplitudes  (i.e.,  small  driving  pressures),  and  the  model 
yielded  nonsensical  results  when  it  was  extended  to  higher  transmitted  pressures, 
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precluding  meaningful  comment  at  this  time  about  what  happens  with  these  stronger 
driving  signals. 

For  a  bubble  in  water  the  line  loss  is  estimated  to  be  ~31  dB  for  a  6  /x  m-diameter 
bubble  and  -34.5  dB  for  a  0.9  /xm-diameter  bubble  at  transmitted  sound  pressure 
amplitudes  of  up  to  0.5  atm. 

Combining  the  predicted  line  loss  with  the  noise  measutements  of  section  VI,  we 
estimate  the  S/N  ratio.  At  a  setting  of  3  V  P-P  on  the  3325 A  signal  generator  the  S/N 
ratio  is  estimated  to  be  6  dB  for  a  6  m  m-diameter  bubble  and  only  -2.5  dB  for  a  0.9  n  m- 
diameter  bvbble.  This  is  uncomfortably  low  There  is  no  further  improvement  in  the 
S/N  ratio  to  be  gained  by  signal  averaging,  because  extensive  signal  averaging  has  been 
used  already  when  obtaining  the  noise  measurements.  However,  it  is  obvious  from  the 
sample  data  from  section  III  that  the  noise  level  can  be  reduced  by  at  least  an  order  of 
magnitude  by  subtracting  a  judiciously  chosen  background  measurement  from  the  signal, 
which  gives  hope  that  the  system  may  be  sensitive  enough  to  detect  individual  bubbles. 

If  the  S/N  ratio  proves  unacceptably  low  in  practice,  then  for  a  remedy  we  may 
consider  that  the  simulation  results  show  that  the  signal  level  from  a  bubble  is  roughly 
proportional  to  the  bubble  diameter.  Therefore,  the  S/N  ratio  could  be  increased, 
without  replacing  the  existing  electronic  components  with  more  linear  components,  by 
investigating  larger  bubbles.  This  approach  would  require  that  we  replace  the  existing 
transducers  with  a  pair  designed  for  a  lower  frequency  range,  which  would  encompass 
the  resonance  frequencies  of  larger  bubbles. 
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IX.  CONCLUSIONS 


The  bubble  detector  lias  been  redesigned  and  provided  with  appropriate  software.  It 
now  has  the  potential  to  be  of  use. 

Simulations  of  bubble  oscillations  have  yielded  important  insights.  These 
simulations,  combined  with  measurements  of  the  system  noise,  suggest  that  it  may  be 
difficult  to  detect  a  small  (<;  6  diameter)  bubble  whose  main  resonance  frequency  is 
within  the  range  over  which  the  transducers  are  sensitive.  A  separate  simulation  of 
bubble  size  dynamics  suggests  that  such  small  bubbles  would  be  difficult  to  stabilize  for 
the  purpose  of  calibration.  We  also  have  learned  via  simulation  that  the  signal  strength 
increases  linearly  with  bubble  diameter  and  that  bubble  stability  improves  with  size. 
Obviously,  it  may  be  advantageous  to  replace  the  present  transducers  with  a  pair  that 
functions  over  a  frequency  range  that  encompasses  the  resonance  frequencies  of  larger 
bubbles,  and  we  are  investigating  that  possibility.  In  the  meantime,  experiments  will  be 
undertaken  to  determine  whether  the  existing  system  is  capable  of  ex  vivo  detection  of 
gas  bubbles  in  dog  spinal  cords  during  and  after  decompression. 
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APPENDIX:  HP  Basic  Program  to  Control  HP  3325 A  Signal  Generator 
and  HP  3561A  Signal  Analyzer 


10  !  THIS  PROGRAM  LISTS  THE  INITIAL  SETTINGS  FOR  THE  3325A  SIGNAL 

20  !  GENERATOR,  THE  HP3561A  SIGNAL  ANALYZER,  AND  THEN  OFFERS  THE 

30  !  OPPORTUNITY  TO  CHANGE  SOME  OF  THESE  DEFAULT  SETTINGS.  IT  ALSO 

40  !  TRIGGERS  THE  EXPERIMENT,  EXTRACTS  THE  DATA  FROM  THE  TIME  BUFFER  IN 

50  !  THE  356 1A,  CONVERTS  THE  DATA  TO  AN  AMPLITUDE- VERSUS-TIME  SPECTRUM 

60  !  (WHICH,  BECAUSE  OF  THE  NATURE  OF  THE  EXPERIMENT  IS  ALSO  AN 

7G  !  AMPLITUDE- VERSUS-FORCING  FREQUENCY  SPECTRUM),  AVERAGES  IT  OVER  A 

80  !  USER-SELECTABLE  NUMBER  OF  REPEAT  MEASUREMENTS,  SUBTRACTS  THE 

90  !  "BACKGROUND"  (OR  BASELINE)  MEASUREMENT,  STORES  THE  DATA 

100  !  FOR  LAYER  USE,  AND  PLOTS  IT  TO  BOTH  THE  SCREEN  AND  A  PRINTER. 

no  ! 

120  OPTION  BASE  1  !  MINIMUM  NUMBER  OF  ELEMENTS  IN  AN  ARRAY  IS  1 

130  PRINTER  IS  CRT  !  PRINTS  TO  SCREEN 

140  RANDOMIZE  !  CHANGE  SEED  OF  RANDOM  NUMBER  GENERATOR 

150  ! 

160  INTEGER  Zero(514)  !  THIS  ARRAY  TO  BE  USED  TO  FILL  THE  BUFFER 

170  MAT  Zcro=  (0)  !  AT  I/O  PATH  ’@Tag’  WITH  ZEROES 

180  !  IN  BETWEEN  TRACE  DUMPS  FROM  THE  SIG  ANALYZER. 

190  ! 

200  INTEGER  Tag _field(350),Pange,Rec_size 
210  INTEGER  Raw_data(30,1024)  BUFFER 
220  INTEGER  Raw“raw_data(1024)  BUFFER 

230  INTEGER.  Rawjracc(401)  !  FOR  READING  TRACE  FROM  BUFFER 

240  INT  EGER  N, Limit, Sweeps, Sweep_back,Ii, 18 

250  REAL  Amplitude(512),Backgrnd(512),Negarea,Posarea 

260  REAL  Trace  data(401)  !  FOR  CONVERTING  BINARY  TRACE  DATA  TO  VOLTAGES 

270  DIM  Name$[‘i5],FilcS[30]  !  FOR  NAMING  THE  OUTPUT  FILES 
280  ! 

290  ASSIGN  @Anz  TO  711  !  OPEN  THE  I/O  PATHS 

300  ASSIGN  ©Siggen  TO  717 

310  ASSIGN  (©Header  TO  BUFFER  [350];FORMAT  OFF 
320  ASSIGN  ©Time  buff  TO  BUFFER  Raw_raw_data(*);FORMAT  OFF 
330  ASSIGN  @Tag  fc  BUFFER  [1028J;FORMAT  OFF 
340  ! 

350  GOSUB  Setdefaults 
360  GOSUB  Display_default 
370  GOSUB  Change_default 
380  GOSUB  Init_siggen 
390  GOSUB  Initanz 
400  ! 

410  ! 

420  BEEP  (INT(RND*5127)+81),1 
430  GOSUB  How_many 
440  GOSLIB  Index 
450  ! 

460  NameS  =  "\HPBLP\DATA\" 

470  FileS  -  Namc$&"STRT"&  VAL$(Istart)&".BCK" 

480  BEEP  (INT(RND*5127)+81),1 


490  ! 

500  !  MEASURING  BACKGROUND 
510  ! 

520  INPUT  "HIT  ’ENTER’  TO  MAKE  A  BACKGROUND  MEASUREMENT  ",DummyS 
530  PRINT "" 

540  PRINT  "NOW  COLLECTING  THE  BACKGROUND  MEASUREMENT" 

550  PRINT  "H 

560  N  =  0  !  INDEX  N  =  0  IF  MAKING  BACKGROUND  MEASUREMENT, 

570  !  N  =  1  IF  COLLECTING  REGULAR  DATA 

580  OUTPUT  @Anz;“NAVG";Sweep  back  !  TRACE  DATA  ARE  TIME  AVERAGED  OVER 

590  “  !  "Sweepback"  NUMBER  OF  SWEEPS 

600  OUTPUT  @Anz;"SCAL"  !  3561A  SELF-CALIBRATION 

610  WAIT  3  !  WAIT  FOR  3561A  SELF-CAL 

620  OUTPUT  @Anz;"STRT  !  START  TIME  AVERAGE 

630  FOR  Ii  =  1  TO  Sweep  back 

640  GOSUB  Triggcr  system 

650  GOSUB  Read_data 

660  WAIT  Sweep  tirae  !  WATT  FOR  SIG  GEN  TO  CYCLE  BACK 
670  NEXT  Ii 

680  GOSUB  Convcrt_and_sum 
690  MAT  Backgrnd-  Amplitude 
700  GOSUB  Storedata 
710  GOSUB  Integrate 
720  GOSUB  Plot  data 

730  DUMP  GRAPHICS  !  PIXEL-BY-PIXEL  DUMP  OF  SCREEN  TO  EXTERNAL  PRINTER 
740  ! 

750  CONTROL  @Tag,3;l  !  MAKE  SURE  BUFFER  POINT  WRITER  IS  ON  1ST  BYTE 
760  OUTPUT  @Tag;Zerc(*)  !  FILL  BUTTER  WITH  514  INTEGER  ZEROES 
770  ! 

780  GOSUB  Read  trace  !  DUMP  TRACE  FROM  SIG  ANALYZER  TO  BUFFER 
790  GOSUB  Scaic_trace 
800  GOSUB  Plottrace 

810  DUMP  GRAPHICS  !  PRINTING  TRACE  BY  DUMPING  SCREEN  TO  PRINTER 
820  ! 

830  !  COLLECTING  DATA 

840  ! 

850  N  =  1  !  INDEX  N  =  0  IF  MAKING  BACKGROUND  MEASUREMENT, 

860  l  N  =  1  IF  COLLECTING  REGULAR  DATA 

870  OUTPUT  @Anz;"NAVG';Sweeps  '  TRACE  IS  TIME-AVERAGED  OVER  "Sweeps" 

880  !  NUMBER  OF  SWEEPS 

890  BEEP  (INT(RND*5127U  81),1 
900  ! 

910  PRINT  "" 

920  INPUT  "HIT  ’ENTER’  TO  START  DATA  COLLECTION",Dummy$ 

930  OUTPUT  @Anz;*STRT  !  START  TIME  AVERAGE 

940  FOR  I8  =  Istart  TO  (1000+Istarl)  '  CAN  DO  UP  TO  1000  COLLECTION  CYCLES 

950  FileJ  =  NameS&"AMP"&  VALS(18)&".DAT" 

960  PRINT  ** 

970  PRINT  "" 

980  BEEP  (INT(RND*512/)+81„l) 

99u  PRINT  "NOW  ON  COLLECTION  CYCLE  NUMBER  ";I8 
1000  PRINT  - 
1010  PRINT "" 

1020  OUTPUT  @Anz;"SCAL"  !  SIG  AN  CALIBRATES  ITSELF 


1030  WAIT  3  !  WAIT  FOR  3561A  SELF-CAL 

1040  OUTPUT  @Ad7;"STRT’  !  START  TIME  AVERAGE 

1050  FOR  Ii-1  TO  Sweep 
1060  GOSUB  Trigger  system 

1070  GOSUB  Read_datay 

1080  WAIT  Sweep  time  !  WAIT  FOR  SIG  GEN  TO  CYCLE  BACK  TO  "STAR T 
1090  NEXT  Ii 

1100  GOSUB  Convert_and_sum 
1110  ! 

1120  MAT  Amplitude  =  Amplitude -Backgrnd  !  SUBTRACTING  THE  "BACKGROUND"" 
1130  !  MEASUREMENT  BEFORE  PROCEEDING 

1140  GOSUB  Store_data 
1150  BEEP  (INT(RND*5127)  +  81),1 
1160  ! 

1170  GOSUB  Integrate 
1180  GOSUB  Plot  data 
1190  DUMP  GRAPHICS 
1200  ! 

1210  CONTROL  ©Tag, 3;  1 
1220  OUTPUT  @Tag;Zero(+) 

1230  ! 

1240  GOSUB  Read_trace 
1250  GOSUB  Scale_trace 
1260  GOSUB  Plot  trace 
1270  DUMP  GRAPHICS 
1280  ! 

1290  PRINT  "CYCLE  COMPLETE;  PAUSING  BEFORE  NEXT  CYCLE" 

1300  PRINT  "" 

1310  PRINT  "" 

1320  PRINT  "" 

1330  WAIT  Waittime 
1340  NEXT  18 
1350  ! 

1360  ASSIGN  ©Path  TO  *  !  CLOSE  THE  I/O  PATHS 

1370  ASSIGN  ©Header  TO  * 

1380  ASSIGN  ©Time  huff  TO  * 

1390  ASSIGN  ©Tag  TO  * 

1400  ! 

1410  STOP 
1420  ! 

1430  Set  defaults:  ! 

1440  !  ~  SET  DEFAULTS  FOR  PARAMETERS 

1450  ! 

1460  !  DEFAULTS  FOR  SIGNAL  GENERATOR 
1470  ! 

1480  Sweep  _start_f~l  !  START  SWEEP  AT  1  MHz 
1490  Strt_sweep_unit$  =  "MH"  !  MEGAHERTZ 
1500  ! 

1510  Sweep  _stop_f=  10  !  STOP  SWEEP  AT  7 

1520  Stop_sweepJJnit$  =  MMH"  !  MEGAHERTZ 
1530  ! 

1540  Signal_amp  -  3  !  AMPLITUDE  IS  3 

1550  Sig_amp_unit$  =  "VO"  1  VOLTS 
1560  ! 


SET  DEFAULTS  FOR  ANALYZER 


1570  ! 

1580  ! 

1590  Start  freq_anz=  1  !  START  OF  FREQUENCY  SPAN  IS  AT  0,8 
1600  Strt  f_anz_unit$ - "KHZ"  !  KILOHERTZ 
1610  ! 

1620  Span  frcqjinz=2  !  SPAN  OF  FREQUENCIES  IS  2 
1630  SpanJ  anz_unit$-"KHZ"  !  KILOHERTZ 
1640  ! 

1650  Sweep_time  =  .4/Span_freq_anz  !  SWEEP  TIME  DEPENDS  ON  SPAN 
1660  ! 

1670  Rccsize  =  1 
1680  ! 

1690  RETURN 
1700  ! 

1710  Display_default:  ! 

1720  !  SHOW  CHANGEABLE  DEFAULTS 

1730  ! 

1740  PRINT  "YOU  ARE  CONTROLLING  THE  HP  3325A  SIG  GEN,  THE  HP  3561A  SIG  AN," 

1750  PRINT  "AND  THE  REST  OF  THE  BUBBLE  DETECTOR  FROM  THIS  PC." 

1760  PRINT"" 

1770  PRINT"" 

1780  PRINT  "SOME  OF  THE  DEFAULT  SETTINGS  ON  THESE  DELICATE  INSTRUMENTS" 

1790  PRINT  "ARE  HARD-CODED  IN  THIS  IMMENSELY  COMPLICATED  PROGRAM,  WHICH" 

1800  PRINT  "SHOULD  BE  TOUCHED  ONLY  BY  HIGHLY  TRAINED,  AUTHORIZED  PERSONNEL." 
1810  PRINT  "HOWEVER,  OTHER  SETTINGS  CAN  BE  SELECTED  EVEN  BY  AN" 

1820  PRINT  "IGNORANT  LAYMAN  LIKE  YOURSELF." 

1830  PRINT  TOR  INSTANCE,  THE  START  AND  STOP  FREQS  ^OR  THE  SWEEP," 

1840  PRINT  "THE  TIME  TO  COMPLETE  A  SWEEP,  AND  THE  AMPLITUDE  OF  THE  SWEEP" 

1850  PRINT  "CAN  BE  CHANGED.  IN  ADDITION,  FOR  TIE  ANALYZER,  THE  START  FREQk 
1860  PRINT  "AND  FREQ  SPAN  CAN  BE  CHANGED." 

1870  PRINT  "NOT  ONLY  THAT,  BUT  YOU  CAN  SELECT  HOW  MANY  SWEEPS  WILL  7;EH 
1880  PRINT  "U ME- AVERAGED  INTO  EACH  CYCLE,  AND  HOW  LONG  THE  SYSTEM- 
1890  PRINT  "PAUSES  BETWEEN  CYCLES.  WHY,  IT’S  LIKE  LIVING  IN  SPACE  !!!N 
1900  PRINT"" 

1910  PRINT"" 

1920  INPUT  "PLEASE  HIT  THE  ENTER  KEY",Dummy$ 

1930  RETURN 

1940  Display_setting:  ! 

1950  PRINT"" 

1960  PRINT  THE  SETTINGS  FOR  THESE  INSTRUMENTS  ARE." 

1970  PRINT"" 

1980  PRINT  "FOR  THE  SIGNAL  GENERATOR" 

1990  PRINT"" 

2000  PRINT  "  1  SWEEP  START  FREQ  =";S weep  start  f;Strt  sv/ecp_unit$ 

2010  PRINT  "  2  SWEEP  STOP  FREQ  =r;Swcep_stop_f;Stop_sweep_unit$ 

2020  PRINT"  3  SWEEP  TIME  =";Sweep_timc;MSECONDSH  ~ 

2030  PRINT"  4  SIGNAL  AMPLITUDE  =";Signal  amp;Sig_ampjinit$ 

2040  PRINT  "" 

2050  PRINT  "FOR  THE  ANALYZER" 

2060  PRINT"" 

2070  PRINT  "  5  START  FREQ  =  ";Siart  freq_anz;Strt_f_anz  unitS" 

2080  PRINT"  6  FREQ  SPAN  =  ";Span  freq_anz;Span_f_an7_umt$ 

2090  PRINT  "" 

2100  RETURN 


21 10  ! 

2120  Change_dcfault:  ! 

2130  !  CHANGE  ANY  OF  THE  DEFAULTS????? 

2140  ! 

2150  GOSUB  Display  setting 
2160  ! 

2170  Firstqucstion:  ! 

2180  ! 

2190  INPUT  "DO  YOU  WISH  TO  CHANGE  ANY  OF  THE  SETTINGS  (Y  OR  N)  ?",Reply$ 

2200  SELECT  ReplyS 
2210  CASE  NVn" 

2220  RETURN 
2230  CASE"Yn/Y 
2240  GOTO  Make_changes 
2250  CASE  ELSE 

2260  PRINT  "PLEASE  RESPOND  WITH  Y  OR  N" 

7270  GOTO  First_question 
2280  END  SELECT 
2290  ! 

2300  Make_changcs:  ! 

2310  !  MAKE  CHANGES  IN  SOME  PARAMETERS 

2320  ! 

2330  INPUT  "PLEASE  INDICATE,  (BY  NUMBER)  WHICH  PARAMETER  SHOULD  BE  CHANGED" 
1340  SELECT 
2350  ! 

2360  CASE  1 

2370  ! 

2380  INPUT  "PLEASE  ENTER  THE  NEW  SWEEP  START  FREQ", Sweep  start  f 
'2390  PRINT  "WHAT  ARE  THE  UNITS  (CHOOSE  FROM  THE  FOLLOWING/1 
2400  PRINT "" 

2410  PRINT "  HZ  HERTZ" 

2420  PRINT "  KH  KILOHERTZ- 

2430  INPUT  ■  MH  MEGAHERTZ", Strt_sweep_unit$ 

2440  GOTO  Changcdefault 
2450  ! 

2460  CASE  2 
2470  ! 

2480  INPUT  "PLEASE  ENTER  THE  NEW  SWEEP  STOP  FREQ",Sweep_stop  f 
2490  PRINT  "WHAT  ARE  THE  UNITS  (CHOOSE  FROM  THE  FOLLOWING)" 

2500  PRINT  - 

2510  PRINT "  HZ  HERTZ- 

2520  PRINT "  KH  KILOHERTZ- 

2530  INPUT "  MH  MEGAHERTZ-,Stop  sweep  unit$ 

2540  GOTO  ChangedetauU 
2550  ! 

2560  CASE  3 
2570  ! 

2580  INPUT  "PLEASE  ENTER  THE  NEW  SWEEP  TIME,  IN  SECONDS", Sweep  time 

2590  !  IF  CHANGE  SWEEP  TIME,  CHANGE  SPAN,  TOO 

2600  Spanfreq^anz  =  .4/Sweep_time 

26 10  Span_f_anz_unit$ = "KHZ- 

2620  GOTO  Changcdcfault 

2630  ! 

2640  CASE  4 


2650  ! 

2660  INPUT  -PLEASE  ENTER  THE  NL  W  SIGNAL  AMPLITUDE -,S  ignalamp 


2670 

2680 

PRINT  "WHAT  ARE  THE  UNITS  (CHOGSE  FROM  THE  FOLLOWING) 
PRINT  " 

2690 

PRINT" 

VO 

VOLTS" 

2700 

PRINT 

MV 

MILLIVOLTS' 

2710 

PRINT  " 

VR 

VOLTS  RMS" 

2720 

INPUT" 

MR 

MV  RMS",Sig  amp  unitS 

2730  GOTO  Change  default 
2740  ! 

2750  CASE  5 
2760  ! 

2770  INPUT  "PLEASE  ENTER  THE  NEW  START  FREQ  FOR  THE  ANALV2ER-,S:art_frcq_anz 
2780  PRINT  "WHAT  ARE  THE  UNITS  (CHOOSE  FROM  THE  FOLLOWING)" 

2790  PRINT "" 

2800  PRINT H  HZ  HERTZ- 

2810  PRINT H  KHZ  KILOHERTZ- 

2820  INPUT "  MHZ  MEGAHERTZ", Strt  f  anz  unitS 

2830  GOTO  Changedefault 

2840  ! 

2850  CASE  6 
2860  ! 

2870  INPUT  "PLEASE  ENTER  THE  NEW  SPAN  FREQ  FOR  THE  ANALYZER", SpanfreqLanz 
2880  PRINT  "WHAT  ARE  THE  UNITS  (CHOOSE  FROM  THE  FOLLOWING)" 

2890  PRINT  "" 

2900  PRINT  ■  HZ  HERTZ" 

2910  PRINT "  KHZ  KILOHERTZ" 

2920  INPUT "  MHZ  MEGAHERTZ", Span  f  anz  unitS 

2930  !  IF  CHANGE  SPAN,  CHANGE  SWEEP  TIME,  TOO: 

2940  IF  Span_f_  anzunitS  -  "HZ"  THEN  Sweep_time  =  400/Span_freq_anz 
2950  IF  Span  f_anz  unitS  = "KHZ"  THEN  Sweep  time-, 4 /Span  frcqjmz 
2960  IF  Span_f_anz_unitS  =  HMHZ'  THEN  Sweep  time-. 0004/Span  freq_anz 
2970  GOTO  Change^deiault 
2980  ! 

2990  CASE  ELSE 
3000  ! 

3010  PRINT  "PLEASE  CHOOSE  AN  INT  EGER  FROM  1  TO  6H 
3020  GOTO  Change  default 
3030  ! 

3040  END  SELECT 
3050  RETURN 
3060  l 

3070  Init  siggen:  ! 

3080  !  ~  INITIALIZE  THE  SIGNAL  GENERATOR 

3090  ! 

3100  OUTPUT  @Siggen;MST";Swccpj4art_f;5trL_swcep_unit$ 

3110  OUTPUT  @Siggcn;"SP";Swcep  stop_f;Stop_sweep_unit$ 

3120  OUTPUT  @Siggen;"Tr;Swecpjime?SEH 

3130  OUTPUT  @Siggen;"Alvr;Sigual_amp;Sig  amp  unitS 

3140  OUTPUT  @Siggen;"RF2" 

3150  ! 

3160  RETURN 
3170  ! 

3180  lnit  anz;  ! 


3190  !  INITIALIZE  THE  ANALYZER 

3200  !  SEE  HP  MANUAL  FOR  THE  MEANINGS  OF  THE  FOLLOWING  HP-IB  MNEMONICS 
3210  OUTPUT  ©Auz;TACMM  !  TIME  CAFTURE  MODEL 
3220  OUTPUT  @Anz;"SF;Startjrcq_anz;Strt_f_anz_unk$ 

3230  OUTPUT  @Anz;"SP";Span  fieq_anz;Span_f_anz_imit5 

3240  OUTPUT  @Anz;"AVTD  ‘  !  USING  TIME  AVERAGING  ON  THE  3561A 

3250  OUTPUT  @Anz;"FLA T  !  CHOOSE  FLATTOP  WINDOWING 

3260  OUTPUT  @Anz;"SLTA"  !  DEFINE  T  RACE  ’A* 

3270  OUTPUT  @Anz;"VSFS.3MVRMM  !  SET  VERTICAL  SCALE  FOR  TRACE 
3280  !  DISPLAY  ON  FRONT  PANEL  Of  3561A 

3290  OUTPUT  @Anz;"VSLI"  !  LINEAR  VERTICAL  SCALE 

3300  OUTPUT  @Anz;"MAG"  !  CHOOSE  MAGNITUDE  DISPLAYED 

3310  OUTPUT  @Anz;"SLTR"  !  DEFINE  TRACE  *B* 

3320  OUTPUT  @Anz;TIRE"  '  REAL  TIME  DISPLAY 

3330  OUTPUT  ©AnzfTBNRIREC'  !  ONE  TIME  RECORD  AT  A  TIME  IN  THE  BUFFER 
3340  OUTPUT  @Anz;TRGRH 
3350  OUTPUT  ©AnzfEXT* 

3360  OUTPUT  @Anz;"SLOPNEG"  !  FXTERNAL  TRIGGER;  TRIGGERED  BY  NEGATIVE  SLOPE 
3370  ! 

3380  PRINT  "WAITING  WHILE  .ANALYZER  IS  BF^NG  INITIALIZED" 

3390  WAIT  5 
3400  ! 

3410  !  DO  ONE  AUTORANGE  ON  THE  3561A  SIGNAL  ANALYZER: 

3420  OUTPUT  ©Siggcn;"SC"  !  START  CONTINUOUS  RUN  ON  3325A 

3430  WAIT  .5 

3440  OUTPUT  @Anz;"SARG"  !  DO  ONE  AUTORANGE  ONLY 

3450  PRINT"" 

3460  PRINT  "WAITING  FOR  SIGNAL  ANALYZER  TO  AUTORANGE" 

3470  PRINT"" 

3480  ! 

3490  WAIT  8  !  WAIT  FOR  COMPLETION  OF  AUTORANGING 

3500  OUTPUT  ©AnzfSCAF"  !  INTO  TIME  CAPTURE  MODE 
3510  WAIT  Sweep  time +  .5 

3520  OUTPUT  @Anz;HDTBBH  !  DUMPING  TIME  BUFFER  IN  BINARY  TO  PC 
3530  TRANSFER  @Anz  TO  @Hcadcr;COUNT  350, WAIT 
3540  CONTROL  @Header>5;35  !  GET  RANGE  FROM  HEADER 
3550  ENTER  @Header;Raugc 

3560  Factor  -  4/3*  10A((Range  +  4.812) /20) /32768  !  FOR  CONVERTING  TIME  BUFFER  DATA 
3570  !  FROM  BINARY  TO  VOLTAGE 

3580  ! 

3590  OUTPUT  ©  AnzfACAL  OFF  !  DISABLES  AUTOCALIBRATION  OF  3561A 
3600  OU 1  PUT  @Anz;"ARNG  OFF  !  DISABLES  AUTORANGING  OF  356 1A 

3610  OUTPUT  @Siggert;"SSH  !  RESETS  3325A  SIG  GEN 
3620  RETURN 
3630  ! 

3640  How  niany:  !  HOW  MANY  SWEEPS  TO  PERFORM  PER  CYCLE 

3650  !  AND  HOW  LONG  BETWEEN  CYCLES 

3660  ! 

3 '70  PRINT"" 

3689  PRINT"" 

3690  INPUT  "HOW  MANY  SWEEPS  PER  CYCLE  ... 

...DURING  REGULAR  DATA  COLLECTION  (1  TO  20)  ?", Sweeps 
3700  IF  Sweeps>20  OR  Swecpsd  THEN  How_many 
3710  ! 


3720  How  many  back:  !  HOW  MANY  SWEEPS  DURING  BKGD  MEASUREMENT 
3730  INPUT  "HOW  MANY  SWEEPS  PER  CYCLE  DURING  ... 

...BACKGROUND  MEASUREMENT  (1  TO  30)  ?"tSwecp  back 
3740  IF  Sweep  back>30  OR  Sweep  back<l  THEN  How  many  back 
3750  ! 

3760  Howlong:  ! 

3770  INPUr  "HOW  LONG  (IN  SECONDS)  BETWEEN  CYCLES  ?",Waittimc 
3780  IF  Waittime  <10  THEN 

3790  PRINT  TOO  SHORT  AN  INTERVAL,  CHOOSE  A  TIME  OF  AT  LEAST  10  SEC- 
3800  GOTO  How  long 
3810  END  IF 
3820  ! 

3830  PRINT  " 

3840  RETURN 
3850  ! 

3860  Index:  » 

3870  INPUT  THE  DATA  FILES  WILL  BE  NUMBERED,  STARTING  WITH  NUMBER  ~",Istart) 
3880  PRINT"" 

3890  RETURN 
3900  ! 

3910  Trigger  system:  !  BEGIN  THE  DATA  COLLECTION 
3920  ! 

3930  OUTPUT  @Anz;"SCAP"  !  ANALYZER  IN  START  CAPTURE  MODE 
3940  WAIT  Sweep  time  +  .5 

3950  OUTPUT  @Siggcn;"SSSS"  !  THE  FIRST  "S5"  RESETS  THE  SIGNAL  GEN 

3960  !  AND  THE  SECOND  STARTS  THE  SWEEP 

3970  WAIT  Sweep  time  !  WAIT  FOR  SWEEP  TO  FINISH  BEFORE 

3980  !  TRYING  TO  READ  TIME  BUFFER 

3990  RETURN 

4000  ! 

4010  Read  data:  i  GET  BUFFERED  DATA  FROM  THE  ANALYZER 

4020  ! 

4030  IF  N  =  0  THEN  Limit  =  Sweepback  !  N  =  0  MEANS  BKGD  BEING  MEASURED 

4040  IF  N -1  THEN  Limit -Sweeps  !  N-l  MEANS  REGULAR  DATA  BEING  COLLECTED? 

4050  PRINT"  ON  SWEEP  OF  ";Limit 

40(H)  PRINT"" 

4070  ! 

4080  OUTPUT  <aAnz;"DTBBM 
4090  New_read:  !" 

4100  CONTROL  ©Header, 3;1 
4110  CONT  ROL  ©Time  buff,3;l 

4120  TRANSFER  @Anz  TO  <o>Header;COUNT  350, WAIT 
4130  TRANSFER  @Anz  TO  @Time_buff, COUNT  2048, WAIT 
4140  MAT  Raw_data(li,*)  =  Rawrawdata 
4150  RETURN 
4160  ! 

4170  ! 

4180  Convert  and  sum:  !  GET  AN  AVERAGE  FOR  THE  DATA  OVER  SEVERAL  SWEEPS 
4190  ! 

4200  PRINT"" 

4210  IF  N  =  0  THEN  !  N  =  0  MEANS  BKGD  BEING  MEASURED 

4220  Limit  =  Sweepback 

4Z*0  PRINT  "NOW  AVERAGING  OVER  THE  ";Limil;"  SWEEPS" 

4240  END  IF 


4250  IF  N  =  1  THEN  !  N  =  1  MEANS  REGULAR  DATA  BEING  COLLECTED 

4260  Limit  =  Sweeps 

4270  PRINT  "NOW  AVERAGING  OVER  THE  SWEEPS  AND  SUBTRACTING  BACKGROUND- 

4280  END  IF 
4290  ! 

4300  FOR  J-l  TO  512 
4310  Amplitudc(J)  =  0 
4320  FOR  Ii=l  TO  Limit 

4330  !  CONVERT  TO  AMPLITUDE  AND  DO  TIME  AVERAGING  SIMULTANEOUSLY 

4340  Amplitudc(J)  =  SQRT(Raw_data(Ii,2*J-l)A24  Raw_data(Ii,2*J)^2)+Amplitudc(J) 

4350  NEXT  Ii 
4360  NEXT  J 

4370  MAT  Amplitude=  Amplitudc*(Factor/Limit)) 

4380  RETURN 
4390  ! 

4400  ! 

4410  Store  data:  !  STORE  AMPLITUDE  SPECTRUM  AFTER  SUBTRACTING  BACKGROUND 

4420  ! 

4430  CREATE  BDAT  File$, 5 12,8 
4440  ASSIGN  ©Path  TO  FilcS 
4450  OUTPUT  @Path;Amp!itudcO 
4460  ! 

4470  RETURN 
4480  ! 

4490  Integrate:  !  FIND  AREA  UNDER  VOLTAGE  AMPLITUDE  SPECTRUM; 

4500  !  SEPARATE  POSITIVE  AND  NEGATIVE  CONTRIBUTIONS  TO  TOTAL  AREA 

4510  Arca  =  0 

4520  Fosarea=0 

4530  Negarea  =  0 

4540  FOR  Irt  =  1  TO  512 

4550  Area  =  (Amplitude(lnt)  +  Ampliludc(lnt  + 1))  +  Area 

4560  IF  Amplitude(Int)>  =0  AND  Amplitude(Inl  +  l)>  -0THEN 

4570  Posarca  -  Posarea  +  Amplitudc(Inl)  +  Amplitude(lnt  + 1) 

4580  ELSE 

4590  IF  Amplitude(Int)  <  =  0  AND  Amplitude(Int  + 1) <  =  0  THEN 
4600  Ncgarea  =  Negarea  + Amplitude(lnt)+Aniplitude(Ini  +- 1) 

4610  ELSE 

4620  IF  Amplitudc(lnt)>0  AND  Amplitude(Int+  1)<0  THEN 

4630  Posarca  -  Amplitude(lnt)*AmpUtude(Int)/(Amplitude(Int)-Amplitude(lnt  + 1))  + Posarca 

4640  Negarca  =  Amplitudc(lnt  +  l)*(l  +  AmpHtude(Int)/(Amplitude(Int  +  l)-Amplitude(Int)))4*Ncgarea 

4650  ELSE 

4660  IF  Amplitude(lnt)<0  AND  Amplitude(lnt  +  1)>0  THEN 

4670  Posarca  =  Amplitudc(Int  + 1)*(1  +  Amplitude(lnt)/(Ainplitudo(lm  +  l)-Amplitude(Int)))  + Posarca 

46S0  Negarca  =  Amplitude(Int)  *  Amplitude(Int)/(Amplitude(lnt)-Amplitude(Int  + 1))  +  Negarea 

4690  END  IF 

4700  END  IF 

4710  END  IF 

4720  END  IF 

4730  NEXT  lnt 

4740  ! 

4750  IF  Strt_swecp  unit$  =  "HZ,<  THEN  Slartf=S\veep_starM/l.E  +  6 
4760  IF  Strt_swecp_unit$  =  "KIT  THEN  Startf=Sweep_starijyi.E  +  3 
4770  IF  St^t_sweep_unit$  =  HMH,,  THEN  Startf=Sweep_slart_f 
4780  IF  Stop_sweep_unit$="HZ’'  THEN  Stopf=Sweep_stopjyi.E+6 


4790  IF  Stop_swccp_unilS  'KH  TIll.N  Stop!  -  1 1*3 

4800  IF  Stop  swecp_unitS  "Mil '  TIILN  Stop!  »  S»ccp_tlof_f 
4810  ! 

4820  Area  =  Arca*.5"(Slo|*f  Siarif)"  1.1.  *  o  512  ’  131NMRT3  TU  il  MM; 

4820  Posarca  -lWrea'^tStopI  StarilVl  L  ♦<>  512  *  CONS LJi  In  TO  »S  Vt«. 

4840  Ncgarca -- NcgjfcaV5*tStof*l-M4rtf)*l  I  ***'512  '  CONS  tK TS  70  *\  MIL 

4850  ! 

46o0  RETURN 
4870  ! 

4880  Pioldata  ’  DISPLAY  ll«.  DATA  ON  sCRLi-N 
4VAI  ! 

4>X*J  CLLA1'  M5RLLN  1  Al.riiA  IHVLA3  KfcTAJNLL)  A.',  OR. AFHKb  D4SPLA3  tS  UALLFL'  l  P 
•OH)  (k.’lJ  < 

4-)20  GISII 
4-0*1  ’ 

4>M0  !  Dl.TLRMiM  MOV*  IAK  Till  S  AAl>  MKK  LD  LXTLNI 

4*)50  RLA1  Ab&4ropx552> 

4‘XjO  MAT  Absamp  ■  A  B8*  Amplitude;  '  OINSLKI  Asptaaic  MATRIX  Tv)  MAuNTH  l-L> 

4oai  Man  •  MAXt  Ahsampt  * »)  *  IH.TIR.MiM  VtAGMU'IH  Of  LaKv.InTS  UlU 

4>>80  Ylim *  ItW’llNT iMxu*  1 .1  -4,-1)  '  MAX  SaUT  ON  >  AXIS,  l  LOT  IS  MIVROSOUS 
4i»i  ! 

5000  IF  N  =  0  THEN  !  '  ABLl  Till.  PLOT  AND  THL  AXl-8 

5010  most:  28,95 

5020  LABEL  'BACKGROUND  <&  "nTIMES.TIMLDATEWD  "aDATUJ  IMLDAT  L) 

50.30  ELSE  !  N  --  i,  SIGNAL  BEING  MEASURED 

5040  MOVL  50,05 

5050  label  -signal  c <1  -atimes<timedatl) 

5060  END  IF 
5070  ! 

5080  MOVE  9,75 
5090  LabelS  =  "Voltage* 

5100  FOR  1  =  1  TO  7 
5110  LABEL  Label$[I,I] 

5120  NEXT  I 
5130  ! 

5140  MOVE  13,80 
5150  LabelS  =  "Amplitude- 
5160  FOR  1  =  1  TO  9 
5170  LABEL  Labcl$[l,I) 

5180  NEXT  I 
5190  ! 

5200  MOVE  41,17 

5210  LABEL  "Forcing  Frequency  (MH/.)" 

5220  MOVE  18,89 

5230  LABEL  VAU(Ylim)&“  uV" 

5240  MOVE  22,58 
5250  LABEL  "0  uV" 

5260  MOVE  15,29 

5270  LABEL.  VALS(-Ylim)&”  uV" 

52SO  MOVE  36,23 
5290  ! 

5300  IFStrt  sweep  unitS  =  "HZ"  THEN  LABEL  VAU(SwecpjstartJ/l.E  +  6) 

5310  IF  Strfsweep”unit$  =  ’KH"  THEN  LABEL  VAL$(Swcepj;tart J/1E  +  3) 

5320  IF  Strt  sweep  unitS  =  "Ml  1"  THEN  LABEL  VALS(Sweep_start_f) 


5330  MOVE  106,23 

5340  IF  Stop_sweep_unit$  =  "H7,”  THEN  LABEL  VAL$(Swcep_stop_f/l.E  +  6) 

5350  IF  Slop_swccp_unit$  =  "KlI"  THEN  LABEL  VALS(Swccp_slop  f/l.E  +  3) 

5360  IF  Slop_oWccji_unitS  =  "MH"  THEN  LABEL  VALS(Swccp_stop_J) 

5370  ! 

5380  MOVE  6,9 

53‘X)  LABEL  USING  ”17A,S6D.2D,15A";"POSlTIVE  AREA  IS  ".Posarca,"  roicrovolls-MHz" 
5400  MOVE  6,4 

5410  LABEL  USING  ‘17A,S6D.2D,15A';”NEGATIVE  ARE\  IS  \Ncgarca,"  miciovolts-MHz" 
5420  ! 

5430  VIEWPORT  38,108,30,90  !  POSITION  THE  GRAPH  ON  THE  CRT 

5440  WINDOW  0, Sweep  timc,(-Ylini/l.E  +  6),(Ylini/l.E  + 6)  !  SIZE  THE  AXES 

5450  AXES  Sweep_time/10,(Ylim/  1.E  +  6)/L0,0,0  !  ADD  TICK  MARKS 

54<^l  ! 


5470  PEN  1 

5480  Dx  -  (Sweep  jime/(Rec_si/e))/5l2 
54‘A)  X=0 

5500  MOM  OAmplitude(l) 

5510  FOR  !  »  1  TO  512 

5520  DRAW  X.AmplituJe(I)  !  PLOT  THE  DATA 

5530  X  X  ■*  Dx 

5540  NEXT  l 

5550  RETURN 

55oO  ! 


5570  Read  trace:  !  READ  THE  MAGNITUDE  TRACE  FROM  THE  3561A, 

5580  !  MAGNITUDE  HAS  BEEN  TIME-AVERAGED  ON  THE  3561A. 

5590  OlTPin  (Vi  Anz;‘SLTA"  !  ACTIVATE  T  RACE  A,  THE  MAGNITUDE  TRACE 
5<i00  OUTPUT  («J Aaz;‘DSTB”  !  DUMP  ACTIVE  TRACE  AND  HFADER  DATA 

5<>H)  TRANSFER  (VFAnz  TO  @Tag;END,WAlT  !  TRANSFER  TO  BUFFER  ON  PC 
5620  CONTROL  (<t  Tag,5;5  !  POSITION  THE  BUFFER  POINT  READER  TO  BYTE  5 

5630  ENTER  @Tag,RawJracc(‘)  !  READ  THE  DATA  FROM  THE  TRACE  INTO  AN  ARRAY 


5640  ! 

5650  RETURN 
5660  ! 


5670  Scale  trace:  !  CONVERT  BINARY  DATA  TO  VOLTAGE 

5680  FOP"  It  =  1  TO  401 


5690  Traccdata(It)  =  10.~(Raw  trace(It)*.005/20.)  !  CONVERTS  TO  VOLTAGE 

5700  NEXT  It 
5710  RETllRN 
5720  ! 


5730  Plot  trace:  !  DISPLAY  THE  TRACE  FROM  THE  3561A  ON  SCREEN 
5740  ! 

5750  CLEAR  SCREEN  !  ALPHA  DISPLAY  RETAINED  AS  GRAPHICS  DISPLAY  IS  CALLED  UP 
5760  GCLEAR 
5770  GINIT 
5780  ! 

5790  Maxi  =  MAX(Trace_dala(")) 

5800  Ymax=100*(lNT(Maxi*l.E  +  4)  +  l)  !  MAX  VALUE  OF  Y;  PLOTTING  IN  MICROVOLTS 
5810  !  DETERMINE  HOW  FAR  THE  Y-AXIS  SHOULD  ^XTEND: 

5820  Ylim  =  Ymax 
5830  ! 

5840  IF  N=0  THEN  !  LABEL  THE  PLOT  AND  THE  AXES 

5850  MOVE  47,95 

5860  LABEL  "BACKGROUND  TRACE" 


5870  ELSE  !  N  =  1;  SIGNAL  BEING  MEASURED 

880  MOVE  50,94 

890  I  ABEL  “SIGNAL  TRACE- 

5900  END  IF 

5910  i 

597-0  MOVE  9,75 
5930  Labcl$  =  "Voll age" 

5940  FOR  1  =  1  TO  7 
5950  LABEL  Labcl$[U] 

5960  NEXT  I 
5970  i 

5980  MOVE  13,80 
5990  Labels  =  "Amplitude" 

6000  FOR  1  =  1  TO  9 
6010  LABEL  Ubcl$[I,I] 

6020  NEXT  I 
6030  ! 

6040  MOVE  42,17 

6050  LABEL  "Offset  Frequency  (kHz)" 

6060  MOVE  18,89 

6070  LABEL  VALS(Ymax)&"  uV" 

6080  MOVE  22,29 
6090  LABEL  "0  uV" 

6100  ! 

6110  IF  Strt_f_anz_  unitS  =  "HZ"  THEN  Stait_freq_an7  =  Start_frcq_anz/l.E  +  3 
6120  IF  Slrt  f  an?  unitS  =  "MHZ"  THEN  Start  freq_anz  =  Start  Jreq_anz*l  .E  +  3 
6130  IF  Span  "f  anz_unit$ = "HZ"  THEN  Span  freq  anz  =  SpanJVcq_anz/LE  +  3 
6140  IF  Span  f  anz_unit$  =  "MlIZ"  THEN  Span_frcq_anz  =  $pan_frcq^anz*l.E  +  3 
6150  1 

6160  MOVE  36,23 

6l7ii  LABEL  VALS (Star tfreqjin 7.) 

6180  MOVE  106,73 

6190  LABEL  VALS(Start_frcq_anz  +  Span  freq^.nz) 

6200  ! 

6210  VIEWPORT  38,108,30,90  !  POSITION  THE  GRAPH  ON  THE  CRT 

6220  WINDOW  0,SpanJreqjmz,0,(YUm/l.E-6)  !  SIZE  THE  AXES 
6230  AXES  Span  frcqjmz/lG,(Ylim/LE  +  6)/10,0,0  !  ADD  TICK  MARKS 

6240  ! 

6250  PEN  1 

6260  Dx  =  Span_frcq_anz/400 
6270  X  =  0 

6280  MOVE  0,Tracc_data(l) 

6290  FOR  It  =  l  TO  401 

6300  DRAW  X,Tracc_data(It)  !  PLOT  THE  DATA 

6310  X  =  X  +  Dx 

6320  NEXT  It 

6330  RETURN 

6340  ! 

6350  END 


